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I.  INTRODUCTION 


Secondary  muzzle  flash  results  from  the  reignition  of  hot,  fuel-rich  gun 
muzzle  exhaust  gases  when  they  mix  with  air  after  the  gun  projectile  is 
launched.  Secondary  muzzle  flash  has  several  deleterious  effects,  so  that 
there  have  been  continuing  efforts  to  learn  to  model,  predict,  and  suppress 
it.  The  Yousefian  flash  prediction  model,  which  includes  the  Muzzle  Exhaust 
Flow  Field  (MEFF)  program^  is  the  only  operational  flash  prediction  model 
that  takes  detailed  chemistry  into  account.  If  one  wants  to  predict  the 
effect  of  a  new  suppressant  combination  or  of  a  new  propellant  composition, 
the  Yousefian  model  is  the  only  game  in  town. 

The  name  MEFF  is  used  by  its  author  to  describe  both  the  front-end 
program  which  models  the  gun  muzzle  exhaust  flow  field  and  the  overall 
modeling  procedure,  which  includes  the  Low-Altitude  Plume  Prediction  (LAPP) 
model.^  LAPP  contains  the  detailed  chemical  modeling;  and  MEFF  was  written  to 
produce  rational  input  parameters  for  LAPP  in  a  fashion  that  changes 
appropriately  with  changes  in  gun  parameters.  MEFF  is  a  ”gun  input”  to  LAPP, 
if  you  would.  In  this  report,  I  shall  attempt  to  confine  my  use  of  the  name 
MEFF  to  descriptions  of  the  muzzle  exhaust  flow  field  program,  with  the  term 
"Yousefian  model"  used  to  describe  MEFF,  LAPP,  and  their  associated  programs. 

The  physics  described  by  MEFF  is  well  documented,  and  the  reader  is 
referred  to  Reference  1  for  questions  concerning  the  reasons  MEFF  is  written 
the  way  it  is.  MEFF  starts  with  the  equations  derived  by  Corner;^  one 
resulting  limitation  of  the  code  is  that  it  is  limited  to  cases  for  which  ttie 
charge  weight  is  significantly  less  than  the  projectile  weight.  Thus,  MEFF 
cannot  be  used  to  model  most  high-velocity  guns,  such  as  those  on  tanks. 

There  are  several  discrete  steps  involved  in  using  the  Yousefian  model, 
and  not  all  of  them  are  obvious.  This  report  is  intended  to  provide  a  well- 
described  path  for  intelligent  MEFF/LAPP  utilization;  one  should  be  able  to 
get  the  desired  results  correctly  and  quickly  by  following  this  guide.  As  an 
aid  to  understanding,  a  sample  calculation  is  followed  from  beginning  to  end, 
and  all  the  input  and  output  are  discussed  thoroughly. 

Four  programs  are  essential  to  the  Yousefian  model.  First,  one  needs  a 
thermodynamics  code;  BLAKE^  is  used  throughout  this  report,  but  NASA-LEWIS 


1.  V.  Yousefian,  "Muzzle  Flash  Onset,"  ARI-RR-236,  Aerodyne  Research,  Inc., 
Billerica,  MA,  November  1980.  Also  available  as  ARBRL-CR-00477,  USA  ARRADCOM, 
Ballistic  Research  Laboratory,  Aberdeen  Proving  Ground,  MD,  February  1982  (AD 
B063  573L). 

2.  R.  R.  Mikatarian,  C.  J.  Kau,  and  H.  S.  Pergament,  "A  Fast  Computer  Program 
for  Nonequilibrium  Rocket  Plume  Predictions,"  AFRPL-TR-72-94,  Air  Force  Rocket 
Propulsion 'Laboratory,  Edwards  Air  Force  Base,  CA,  August  1972. 

3.  J.  Corner,  Theory  of  Interior  Ballistics  of  Guns,  John  Wiley  &  Sons,  New 
York  (1950). 

4.  E.  Freedman,  "BLAKE  -  A  Thermodynamic  Code  Based  on  TIGER:  User's  Guide 
and  Manual,"  ARBRL-TR-02411 ,  USA  ARRADCOM,  Ballistic  Research  Laboratory, 
Aberdeen  Proving  Ground,  MD,  July  1982  (AD  A121  259). 
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could  have  been  used.  This  report  contains  enough  details  for  an  experienced 
BLAKE  user  to  get  the  desired  results;  a  novice  will  probably  have  to  get  help 
using  BLAKE.  The  thermodynamics  code  is  used  twice,  for  two  separate 
functions.  Second,  an  interior  ballistics  code  is  used;  1  have  here  used 
IBHVG,  a  modern  version  of  the  Baer-Frankle  lumped-parameter  model. ^  This 
report  contains  enough  details  for  an  occasional  IBHVG  user  to  get  the  desired 
results;  a  novice  will  have  to  get  help  running  IBHVG.  MEFF  is  needed  as  the 
third  major  step,  and  the  information  in  this  report  is  intended  to  be 
sufficient  to  run  MEFF.  Finally,  LAPP^  is  needed;  and  again,  the  information 
in  this  report  should  be  all  a  user  needs  to  use  LAPP  for  this  application. 


For  ease  in  doing  MEFF  calculations,  I  wrote  two  diort  linking  programs, 
MTOB  and  CONCEN.  They  enabled  automation  of  MEFF  calculations  to  the  maximum 
extent  practicable. 


II.  PRELIMINARY  IHERMOCHEMICAL  CALCULATION 

Required  Results  of  the  Calculation:  First  of  all,  one  must  input 
appropriate  propellant  data  into  a  thermodynamics  code.  The  needed  output  are 
the  quantities  required  by  the  interior  ballistic  code  and  by  MEFF  to  follow. 

Here,  IBHVG  will  be  used  for  the  interior  ballistic  code,  so  one  needs: 

Impetus,  Adiabatic  Flame  Temperature,  Gamma,  and  Covolume 
For  MEFF,  one  also  needs: 

Molecular  Weight 

An  Example  of  the  Calculation:  The  example  that  has  been  chosen  for 

this  illustrative  calculation  is  for  a  155-mm  howitzer;  and  the  propellant  is 
the  standard  M30A1  propellant,  which  contains  Mi  (by  weight)  of  K2SO4  flash 
suppressant.  A  listing  of  the  input  job  stream  and  data  for  this  calculation 
is  included  as  Appendix  A. 

Note  the  deliberate  suppression  of  many  condensed  species,  e.g.,  KCO$, 
KSO$,  K$,  etc.,  for  the  calculations.  Since  the  suppressant  is  presumed  to 
operate  in  the  gas  phase,  solid-phase  or  liquid-phase  final  constituents  that 
could  conceivably  tie  up  some  of  the  potassium  were  not  permitted  to  be  formed 
in  these  calculations. 

The  line  which  begins  with  CM2  is  the  listing  of  the  propellant 
constituents  and  the  weight  percentage  of  each  in  the  propellant: 

NC1260  nitrocellulose,  12.60%  27.90%  of  the  total 


NG  nitroglycerin  22.42% 
NQ  nitroguanidine  46.84% 
EC  ethyl  centralite  1.49% 
KS  potassium  sulfate  1 .00% 


5.  P.  G.  Baer  and  J.  M.  Frankie,  "The  Simulation  of  Interior  Ballistic 
Performance  of  Guns  by  Digital  Computer  Program,”  BRL  Report  No.  1183,  U.S. 
Army  Aberdeen  Research  and  Development  Center,  Ballistic  Research 
Laboratories,  Aberdeen  Proving  Ground,  MD,  December  1962  (AD  299  980). 
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ALC  ethyl  alcohol  0.25$ 

C  carbon  (graphite)  0.01? 

In  this  case,  a  GUN  calculation  is  desired,  with  the  standard  loading 
density  of  0.2  g/cc.  '  3 

The  results  of  this  calculation  are  shown  in  detail  in  Appendix  B.  The 
parameter  values  which  carry  over  directly  to  the  interior  ballistic  code  are: 


Impetus 

Flame  Temperature 

Gamma 

Covolume 


356461  ft-lb/lb 
3003  K 
1.2412  - 
28.81  in^ 


The  values  needed  for  MEFF  are: 


Molecular  Weight  23.432 
Covolume  1.041  cc/g 

III.  INTERIOR  BALLISTIC  CALCULATION 

Input  and  Output  for  the  Calculation  Next,  one  must  use  an  interior 
ballistic  code.  The  input  parameters  needed  are: 

Chamber  volume,  length  of  travel,  propellant  mass,  projectile  mass,  and 
barrel  cross  section 

The  output  quantities  muzzle  velocity  and  mean  gas  temperature  at  shot 
ejection  are  needed  by  MEFF,  as  are  several  of  the  gun  parameters. 

An  Example  af  the  Calculation:  Appendix  C  is  an  IBHVG  calculation  for  the 
example  system.  System,  projectile,  and  propellant  parameters  are  nominal 
values  for  this  system.  Note  the  thermodynamic  values  introduced  from  the 
prior  BLAKE  calculation. 

The  results  of  the  interior  ballistic  calculation  needed  for  MEFF  are: 

Muzzle  velocity  2650  ft/s 

Mean  gas  temp  1860  K 

IV.  MEFF  CALCULATION 


At  last  one  comes  to  the  actual  calculation  with  MEFF  itself.  The  input 
needed  are  illustrated  by  the  data  on  the  bottom  of  the  FLASH  job  stream 
included  as  Appendix  D.  The  MEFF  listing  itself,  as  modified  slightly  for 
automated  running,  is  included  as  Appendix  E.  MEFF  requirements  are  as 
follows: 


Datacard  Requirement 

1  V  A  title  card 

2  Muzzle  velocity 

2  Chamber  volune 

2  Travel 


Illustrative  value 

155-MM  HOWITZER  WIIH  M203  CHARGE 

789  m/s  _ 

.01966  m^ 

5.08  m 
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2  Propellant  mass 

2  Projectile  mass 

2  Bore  cross  section 

2  Gamma 

2'  Molecular  Weight 

2  Mean  gas  temperature 

at  shot  ejection 

2  Covolune 

3  No.  iterations  between 

stored  values 
3  Step  size 

3  Init.  condition  step 

away  from  tau=0 

4  Maximun  distance  from 

muzzle  for  calcula¬ 
tion  to  proceed 
4  Print  step 

4  Diffusion  step  size 

passed  to  LAPP 

4  LAPP  output  parameter 

=1 ,  all  output 
=0,  centerline 
temperatures  only 

I  have  never  changed  the  value  of  the  parameters  on  the  third  data  card.  I 
added  tiie  fourth  card  so  it  would  be  easy  to  vary  the  weapon-dependent  maximum 
distance  from  the  muzzle  for  the  calculations,  the  print  step,  and  the 
diffusion  step  size.  The  fourth  parameter  on  the  card  makes  it  easy  to  reduce 
the  voluminous  LAPP  output  to  just  centerline  temperatures,  for  trouble¬ 
shooting.  The  diffusion  step  size  can  be  increased,  and  the  run  time  will  be 
shorter.  When  one  increases  the  diffusion  step  size  too  far,  the  program  will 
gracefully  halt  at  the  moment  of  ignition.  Too  large  a  diffusion  step  size 
will  not  lead  to  improper  results  of  these  calculations. 

The  output  parameters  passed  from  MEFF  for  calculations  to  follow  are 
written  to  TAPE9  at  statement  5011,  and  they  are  as  follows: 

TN,  the  gas  temperature  at  the  normal  shock  when  the  velocity  of  the  muzzle 
gas  flow  becomes  sonic 

TM0,  the  muzzle  gas  temperature  at  the  time  of  shot  ejection 
TB,  the  gas  temperature  at  the  mixing  region  boundary  when  the  velocity  of 
the  muzzle  gas  flow  becomes  sonic 

PM0,  the  muzzle  gas  pressure  at  the  time  of  shot  ejection 
UB,  the  gas  velocity  at  the  mixing  region  boundary  when  the  velocity  of  the 
muzzle  gas  flow  becomes  sonic 
ALPHA1 ,  the  fraction  of  gas  entering  the  reflected  shock 
RB,  the  radius  of  the  mixing  region  boundary  when  the  velocity  of  the  muzzle 
gas  flow  becomes  sonic 

XMAX,  the  maximum  distance  from  the  muzzle  for  the  calculations,  in  meters 
PRNT,  the  print  step,  in  meters 

FDL,  the  diffusion  step  size,  passed  through  for  LAPP 
KEY,  the  LAPP  output  parameter,  passed  through  for  LAPP 


12.23  kg  (total  of  all) 
46.36  kg 
.0192  m^ 

1.243 
23.43 
1861  K 

.001046  m^/kg 


.001 

.01 


50.  (meters.  10  meters  would  be 
appropriate  for  a  mortar 
calculation. ) 

5.  (meters) 

.2 
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V.  MTOB  CALCULATION 


Next  one  prepares  for  thermodynamic  calculations  at  several  different 
places  in  muzzle  exit  gas  flow  space.  MTOB  (Appendix  F)  was  written  to  do 
this  automatically.  It  takes  the  MEFF  output,  combines  it  with  the  details  of 
the  propellant  contained  in  BOIL  (Appendix  G),  and  produces  TAPE8,  which  is  a 
detailed  command  stream  for  BLAKE  and  data  for  programs  CONCEN  and  LAPP  to 
follow. 

VI,  BLAKE  CALCULATION 

Here  one  needs  two  thermodynamic  calculations,  one  to  calculate  the  mole 
fractions  at  the  normal  ^ock,  and  one  to  calculate  the  mole  fractions  at  the 
reflected  shock.  The  first  of  these  is  simply  a  "point”  calculation  at  the 
pressure  and  temperature  of  the  normal  shock.  For  the  second,  one  recalls' 
that  the  propellant  gas  expands  isotropically  from  the  muzzle  to  the  reflected 
shock  region,  so  that  the  mole  fractions  are  the  same  as  those  at  the  muzzle; 
one  thus  does  a  "point”  calculation  at  the  pressure  and  temperature  of  the 
muzzle  gas  as  it  emerges  from  the  gun. 


VII.  CONCEN  CALCULATION 

The  program  CONCEN  (Appendix  H)  reads  the  BLAKE  output  and  automatically 
calc!j].ates  the  mole  fraction  for  each  gaseous  species  at  the  initial  boundary, 
as  shown  in  Reference  1. 

Xi  =  (1-a)Xjj+a  x^,  where 

Xj^  is  the  mole  fraction  at  the  initial  boundary, 

Xjj  is  the  mole  fraction  at  the  normal  stiock, 

Xp  is  the  mole  fraction  at  the  reflected  shock,  and 

a  is  the  fraction  of  the  flow  that  enters  the  reflected  shock. 

Thus,  the  output  of  CONCEN  or  of  a  calculation  by  some  other  means,  is  a  list 
of  the  13  gaseous  species  that  LAPP  will  consider,  and  the  mole  fraction  of 
each,  in  the  exact  order  that  LAPP  expects  to  find  them:  H2O,  CO,  Hp,  Np, 

cop,  .H,  OH,  0,  Op,  K,  KOH,  KOp,  HOp.  These  results  are  passed  to  Lifpp  on 
TAx  E2  • 

VIII.  LASTDA,  THE  LAPP  STANDARD  DATA  SET 

The  file  LASTDA  contains  the  standard  LAPP  input  data,  including 
thermodynamic  information  on  the  gaseous  species  allowed  and  reaction  rate 
data  on  the  reactions  considered.  The  LAPP  report^  documents  the  needed  LAPP 
input  data  in  detail;  here  we  concentrate  on  much-used  or  often-changed  data 
and  on  changes  from,  the  LAPP  report.  I  have  used  numbers  composed  of  all  9*s 
to  "hold  the  place"  of  values  which  will  be  replaced  in  a  subsequent  read,  in 
order  to  minimize  rewriting  LAPP. 
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1 


LAPP 

Input  Source 

Card  of  input  Via  Input 


1 

FLASH 

TAPES  Title  Card 

On 

card  2,  all 

data  entered  in  15 

2 

LASTDA 

Initial  nunber  of  grid  points 

2 

LASTDA 

Nunber  of  species  (24  max);  here  13 

2 

LASTDA 

Viscosity  option  key;  6=Donaldson/Gray 

2 

LASTDA 

Nunber  of  reactions;  here  25;  LAPP  has  been 
modified  to  handle  up  through  49 

2 

LASTDA 

Three  items  specifying  output  options;  all  0 

2 

LASTDA 

Max  computer  time;  never  used;  here  left  at  200 

2 

LASTDA 

Pressure  option;  never  changed;  here  0 

2 

LASTDA 

Number  of  thermo  entries  for  each  sp>ecies;  here 
22  and  never  changed 

3 

LASTDA 

Signal  frequencies;  left  blank  since  no 
attenuation  calculations  are  desired 

On 

card  4,  all 

data  entered  in  E10.3 

4 

LASTDA 

Initial  value  of  X,  the  distance  from  the 
muzzle,  in  meters. 

4 

RASH 

TAPES  Final  value  of  X-in  meters.  •  ' 

4 

aASH 

TAPES  Print  increment  ib  meters  * .  .i  ■  .b 

Lewis  number,  her.e  always  t  _ 

Prandtl' number,  hfere -always  1  -V  -  - 

4 

LASTDA 

4 

LASTDA 

4 

MEFF 

TAPES  Initial  Boundary  Radius,-,  in  meters 

4 

LASTDA 

Factor  used  to  vary  eddy  viscosity,  here 
always  1 

On 

card  5,  all 

data  entered  in  E10.3 

5 

LASTDA 

Minimum  integration  step  size,  here  always 
.IE-10 

5 

FLASH 

TAPES  Diffusion  step  size,  FDL 

5 

LASTDA 

Pressure  coefficients,  here  always  zero 

On  card  6,  all  data  entered  in  E10.3 


6 

LASTDA 

Pressure  at  initial  value  of  X,  in  atmospheres. 

here  always  1 

6 

MEFF 

TAPES 

Temperature  at  initial  boundary  in  K 

- 

6 

LASTDA 

TAPES 

Ambient  air  temperature  in  K,  here  always  294 

Gas  velocity  at  the  initial  boundary,  in  m/s 

6 

MEFF 

6 

LASTDA 

Ambient  air  velocity  in  m/s,  here  always  3.0. 

Should  never  be  set  to  0. 

6 

LASTDA 

'F  .  Here  always  blank,  so  program  calculates 
this  quantity. 

6 

LASTDA 

Kinetics  cut-off  temperature  in  K,  here  always 

7 

200 

8 

LASTDA 

Mole  fractions  of  ttie  allowed  gaseous  species 

9 

in  ttie  ambient  air 

10 


CONCEN  TAPE2 


10 

11 

12 


Mole  fractions  of  tine  allowed  gaseous  species 
at  the  initial  boundary 


Next  in  LASTDA  come  the  toermodynamic  data  on  the  allowed  gaseous 
species,  in  JANAF- table  format.®  These  are  only  changed  when  new  data  are 
published.  Each  species  name  is  entered  in  A6  format,  and  all  data  are  then 
entered  in  E10.3  format.  The  ordering  of  the  species  in  this  table  determines 
the  order  in  which  LAPP  processes  the  species  in  all  of  its  transactions.  It 
even  specifies  which  species  to  associate  with  the  initial  nianber  densities 
which  have  already  been  read  in  from  this  data  set. 


Last  in  LASTDA  are  the  reaction-rate  data  for  the  25  reactions.  The 
reactions  may  be  listed  in  any  order.  The  data  for  the  C-N-O-H  and  K 
reactions  have  been  fairly  thoroughly  used  and  checked.  At  the  time  this 
report  is  written,  they  are  the  best  available  set  of  reactions  and  the  best 
available  reaction  rates  for  those  reactions;  but  they  are  not  represented 
here  as  being  the  final  correct  description  for  the  reactions. 

The  format  for  each  reaction  follows: 


Colunn 

Item 

Format 

1-6 

Species  A 

A6 

7 

8-13 

+  sign 

Species  B  (or  M) 

A6 

14 

15-20 

+  sign  (if  needed) 

Blank  or  M 

21 

22-27 

=  sign 

Species  C 

A6 

28 

+  sign  (if  needed) 

29-34 

Species  D  (or  M) 

A6 

35 

36-41 

+  sign  (if  needed) 

Species  E  (or  M) 

A6 

42-48 

49-50 

Blank 

Reaction  type,  1  to  10  (see  below) 

12 

51 

Rate  coefficient  type,  1  to  7  (see  below) 

11 

52-59 

A,  Pre-exponential  factor,  cm-molecule-sec  units  E8.2 

60-63 

N,  Temperature  exponent 

F4.1 

64-72 

B,  Activation  energy,  cal/mole 

F9.1 

The 

ten  possible  reaction  types  are  these: 

1 

A  +  B 

t  C  +  D 

2 

A  +  B 

+  n  t  C  +  M 

3 

A  +  B 

t  C  +  D  +  E 

4 

A  +  B 

t  c 

5 

A  +  M 

t  C  +  D  +  M 

6 

A  +  B 

C  +  D 

7 

A  +  B 

+  M  ^  C  +  M 

8 

A  +  B 

^  C  +  D  +  E 

6.  D.  R. 

Stull 

and  H.  Prophet,  ”JANAF  Thermochemical  Tables, 

SECOND  EDITION, 

Dow  Chemical  Ccanpany,  Midland,  MI,  July  1970. 
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9  A  +  B  ^  C  • 

10  A  +  M  ->0  +  0  +  M 


The  seven  possible  rate  coefficient  types  are  these: 

1  k  =  A 

2  k  =  AT”^ 

3  k  =  AT“2 

1 

4  k  =  AT“  i 

5  k  =  A  exp  (B/RT) 

6  k  =  Ar^^expCB/RT) 

3 

7  k  =  Ar  2 

X.  FINAL  RESULTS 

Some  of  the  pages  of  the  final  results  of  the  illustrative  calculation 
are  included  as  Appendix  K.  The  total  calculation  took  78  CPU  seconds  on  MFZ, 
all  but  14.5  seconds  of  which  was  LAPP  execution. 

MEFF  results  are  followed  by  BLAKE  calculations  at  the  two  necessary 
combinations  of  temperature  and  pressure.  These  are  followed  by  the  results 
of  CONCEN,  the  combined  mole  fractions  at  the  initial  boundary. 

Finally  come  the  results  of  the  LAPP  calculation.  First  the  input 
parameters  are  echoed.  Then,  at  each  desired  distance  from  the  muzzle,  as  a 
function  of  the  distance  from  the  centerline  of  the  flight  of  the  projectile, 
there  are  gas  temperature,  velocity,  density,  etc.  and  mole  fractions  for  each 
of  the  allowed  gaseous  consitutents.  Included  in  Appendix  K  are  the  prints 
for  the  flow  from  the  muzzle  (X  =  0  meters),  for  X  =  10  meters  (where  one 
notes  a  maximum  gas  temperature  of  1079  K),  for  X  =  15.3  meters,  and  for  X=  42 
m.  The  prints  at  15.3  meters  are  the  most  interesting,  for  they  show  that  the 
maximum  tonperature  has  risen  to  2018  K,  indicating  that  ignition  has  taken 
place.  If  there  had  been  sufficient  suppressant  to  eliminate  the  flash,  one 
would  have  seen  the  temperature  rise  to  1100  K  or  1200  K  and  then  decline 
slowly,  indicating  that  the  mixture  had  cooled  before  ignition  could  take 
place. 

It  was  noted  earlier  that  one  could  set  KEY  =  0  in  the  input  job  stream 
for  MEFF/LAPP  (Appendix  D),  and  that  then  one  would  have  gotten  only  the 
centerline  temperatures  from  LAPP.  Notice  that  even  then,  by  42  meters  the 
centerline  temperature  has  exceeded  2000  K,  so  that  even  with  a  limited  print, 
ignition  is  unmistakably  indicated. 
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APPENDIX  A 


BLA30A1 ,  INPUT  FOR  A  BLAKE  CALCULATION 


GEK, STMFZ , P6 , T1 20 . BLA30A1 
ACCOUNT, XXXXXX. 

ATTACH , TT, BLAKELIBRARY, IDrELI . 

COPY, TT, TAPE?. 

RETURN, TT. 

REWIND, TAPE?. 

ATTACH, B, BLAKE, IDrELI. 

B. 

EXIT. 

*EOR 

TIT,M30A1 

PRL,CON,2 

REJ , N , K2S04 , C , C2 , CH , CH20 , HN03 
REJ,C(S),K2S04$ 

REJ,K0H$,K02$,K202$ 

REJ , H2S, S20 , S02 , K$, K20, K202 
REJ,KC0$,KS0$,K20$,K$ 

REJ,K2C03$ 

REJ,K2S$ 

UNI, ENG 

CM2,NC1260,2?.9,NG,22.42,NQ,46.84,EC,1.49,KS,1., 

ALC,.25,C,.l 

GUN,.1,.1,.5 

QUIT 


4 
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APPENDIX  B 
BLAKE  CALCULATION 
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TOTAL  GAS  IMOLES/KG)  42.7266  42.6770  42.6341 


PROGRAM  BLAKEf  VERSION  205.11  ** 
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APPENDIX  C 
BHVG  CALCULATION 
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GUN  TYPE:  155-MM  HOWITZER 
CHAMBER  VOLUME:  1200.00  CU  IN 
GROOVE  DIAMETER:  6.200  IN 
GROOVE/LAND  RATIO:  1.660 
TWIST:  ONE  TURN  IN  20.0  CALIBERS 
PRESSURE  GRADIENT:  LAGRANGIAN 
PROJECTILE:  M483A1 


BORE  LENGTH:  200.0  IN 
TIME  STEP:  .100  MS 
LAND  DIAMETER:  6.100  IN 
BORE  AREA:  29.828  SQ  IN 
EXPANSION  RATIO:  6.0 
EROSIVE  COEFF:  .0000500 
PROJ  WT:  102.200  LB 


ENGRAVING  &  FRICTIONAL  RESISTANCE  [KPSI]  VS.  TRAVEL  [IN]: 


TRAVEL:  0.00  .40 

RESISTANCE:  .25  3-35 

PROPELLANT  BLK  POWDR 

WEIGHT  [LB]  .315 

IMPETUS  [FT-LB/LB]  96000. 
FLAME  TEMP  [K]  2000. 

ALPHA  0.0000 

BETA  50.000000 

GAMMA  1.250 

COVOL  [CU  IN/LB]  30.000 

DENS  [LB/CU  IN]  .06000 

GRAIN  TYPE  CORD 

GRAIN  LEN  [IN]  .2000 

GRAIN  WIDTH  [IN]  - 

GRAIN  DIAM  [IN]  .1000 

PERF  DIAM  [IN]  - 

GRAIN  THICK  [IN]  - 

INNER  WEB  [IN]  - 

OUTER  WEB  [IN]  - 

IGNITION  CODE  0 


THRESHOLD  VALUE  0.00000 


1.00  1.55  2.05  4.50  200.00 
4.95  3.63  3.25  2.80  1.50 

M30A1  NC  TUBE 

26.150  .500 

356461.  180000. 

3003.  1553. 

.7000  0.0000 

.003950  30.000000 

1.241  1.250 

28.810  30.000 

.05717  .03400 

7-PERF  SHEET 

.9481  28.0000 

-  1.5000 

.4173  - 

.0338  - 

-  .1250 

.0790  - 

.0790  - 

0  0 

0.00000  0.00000 


STDM203,  A  STANDARD  M203  CHARGE  IN  A  155-MM  HOWITZER 


CONDITIONS  AT: 

MAX  PR 

MUZZLE 

PROP  1 
BURNT 

PROP  2 
BURNT 

PROP  3 
BURNT 

TIME  [MS] 

6.65 

13.69 

1.00 

11.20 

2.10 

BR  PRES  [KPSI] 

47.77 

11.48 

3.01 

19.96 

9.04 

MN  PRES  [KPSI] 

46.02 

11.10 

2.91 

19.26 

8.81 

BS  PRES  [KPSI] 

42.52 

10.34 

2.72 

17.88 

8.35 

MEAN  TEMP  [K] 

2622. 

1860. 

2266. 

2088. 

2530. 

TRAVEL  [IN] 

23.9 

200.0 

.0 

124.4 

.4 

VEL  [FPS] 

1104. 

2650. 

9. 

2377. 

48. 

ACCEL  [G«S] 

11615. 

2530. 

644. 

4593. 

1524. 

FR  BRNT  PROP  1 

1.000 

1.000 

1.000 

1.000 

1.000 

FR  BRNT  PROP  2 

.571 

1.000 

.010 

1.000 

.043 

FR  BRNT  PROP  3 

1.000 

1.000 

.502 

1.000 

1.000 
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APPENDIX  D 


FLASH,  THE  JCL  AND  INPUT  DATA  FDR  MEFF,  MTOB,  BLAKE,  CONCEN,  AND  LAPP 


GEK, STMF7 , PI , T300 .  FLASH 
ACCOUNT, XXXXXXX. 

BEGIN , GETMFA, FILE, LF=A, PF=MEFF, UN=GEK. CREATES  TAPE9 
MAP, OFF. 

FTN,I=A,L=0,R=0,T. 

LGO. 

REWIND, TAPE9. 

BEGIN , GETMFA, FILE, LF=E, PFrWTOB,  UN=GEK. 

BEGIN, GETMFA, FILE,LF=TAPE4,  PF=BOIL,UN=GEK.BLAKE  BOILERPLATE 
FTN ,  B=LG03 , 1=  E,  L=0 ,  R=0 ,  T. 

LG03. 

REWIND,  OUTPUT. 

COPY, OUTPUT, TAPES. 

REWIND, TAPES. FILE  OF  BLAKE,  CONCEN,  AND  LAPP  DATA 
ATTACH, TT,BLAKELIBRARY, ID=ELI. 

COPY, TT, TAPE?. 

RETURN, TT. 

REWIND, TAPE?. 

ATTACH,B, BLAKE, ID=ELI. 

B, TAPES. 

REWIND,  OUTPUT. 

REWIND,  TAPE1. 

COPY, OUTPUT, TAPE1. 

REWIND, TAPE1 .  INPUT  FOR  CONCEN 

BEG  IN, GETMFA, FILE, LF=C, PF=CONCEN,UN=GEK. 

FTN,B=LG01 ,I=C,L=0,R=0,T. 

IGOI.  . 

REWIND,  TAPE2. 

BEGIN , GETMFA, FILE, LF=TAPE3 , PF=LASTDA , UN=GEK. 

REWIND, TAPE3. 

BEGIN, GETMFA, FILE, LF=D, PF=LAPP, UNrGEK. 

FTN,B=[G02, I=D, L=0, R=0, T. 

LG02. 

EXIT. 

155-MM  HOWITZER  WITH  M203  CHARGE 

S0? .? , . 01 9664 ,5 .0S, 1 2.23 ,46 .36 ,.019244 , 1 .241 ,23 .432,1 S60. , .  001041 

4. . 001. .01 

50. . 5.0.0.20.1 
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APPENDIX  E 
MEFF  PRCXJRAM  LISTING 
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ooooo  o  ooo  oooooooooooooooooooooooooo 


PROGRAM  MEFF( INPUT, OUTPUT, TAPE5=INPUT, TAPE6=0UTPUT, TAPES ) 


MAIN:  MUZZLE  FLASH 
MUZZLE  PROPERTIES  VS  TIME 
TYPE  DATA  CARDS  USING  NAFELIST  FORMAT 
TITLE  CARD  FIRST 
*PHYSIC* 

V0=MUZZLE  VELOCITY (M/SEC) 

L=BORE  LENGTH(M) 

MP=MASS  OF  PROPELLANT 
W=MASS  OF  PROJECTILE(KG) 

ArCROSS  SECTIONAL  AREA  OF  B0RE(M**2) 

GAM=RATIO  SPECIFIC  HEATS  PROPELLANT  GAS 
MBARrMEAN  MOL  WEIGHT  PROPELLANT 
TFrFLAME  TEMPERATURE (K) 

TA=AVERAGE  BARREL  GAS  TEMPERATURE(K) 

CVOrCHAMBER  VOLUME(M**3) 

ETArCOVOLUME  FOR  VAN  DER  WAAL'S  EQ.  STATE (M**3/KG) 

*MATH* 

MzNUMBER  OF  ITERATIONS  BETWEEN  STORED  VALUES 
DELTAUrSTEP  SIZE 

DLTAU0=INITIAL  CONDITION  STEP  AWAY  FROM  TAU=0 

REAL  MP,L,MBAR,LAMM,M0,MACHNO,MR,MNP2 
COMMON/WORKA/LAMM, GAM, EPS, M0 , B0 , C0 , D0 , KPMAX 
COMMDN/WORKB/THETA ( 500 ) , Z ( 500 ) , HI ( 500 ) , HI P ( 500 ) , IH  ETAP ( 500 ) , 

1  TTAU(500),FEJECT(500) 

COMMON/WORKC/TEFF, ETA, A, L, MP, V0 , ALMP, R 
COMMONAJORKD/PE(50) ,UE(50) ,TE(50) , RHOE(50) ,TIME(50) , 

1  TEMP(50) ,FKP(50) ,CO2(50) ,CO(50) ,H2O(50) ,H2(50) 

DIMENSION  TITLE(20) 

NAMELIST/PHYSIC/V0 , L, MP, W, A, GAM, MBAR , TF, TA, ETA , CVO 
NAMELIST/MATH/M, DELTAU, DLTAU0 
DATA  IREAD,IWRITE/5,6/ 

DATA  IMAX/500/ ,  IZ THET/0/ , PI/3 .1415 92654/ 

READ  IN  DATA 

READ(IREAD,1111)TITLE 
111  FORMAT(20A4) 

WRITE(9,1111)TITLE 
000  FORMAT(1H1,20A4) 

READCIREAD, »)V0 , CVO, L, MP, W, A, GAM, MBAR, TA, ETA 

CALL  NMLST('PHYSIC‘,0,5,V0,L,MP,W,A,GAM,MBAR,TF,TA,ETA,CVO) 

READ( IREAD, * ) M, DELTAU , DLTAU0 

CALL  NMLSTC ' MATH  ',0,5, M, DELTAU, DLTAU0) 

READ  MAX  DIST,  PRINT  STEP,  DIFFUSION  STEP  SIZE,  AND  KEY;  ALL  FOR  LAPP. 
KEYrI  FOR  ALL  LAPP  PRINTS,  =0  FOR  CENTERLINE  ONLY. 

READdREAD,  *)XMAX,  PRNT,  FDL,  KEY 
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c 

CALBER  =  1000.  «  SQRT(  4. W3. 141 593  ) 

C 

IF(TF.EQ.0.0)  GO  TO  1 
WRITE (IWRITE, 1000)  TITLE 

WRITE ( IWRITE ,1010)  V0 , L , MP , W, CALBER , A , GAM , MBAR , TF, ETA , 

1CVO 

1010  FORMAT('0','  MUZZLE  VELOCITY  V0  =»,G10.4,»  M/SECV 
1*  •,»  BORE  LENGTH  L  =',G10.4,’  METERS'/ 

2'  PROPELLANT  MASS  MP  =',G10.4,'  KG'/ 

3'  ','  PROJECTILE  MASS  W  =',G10,4,'  KG'/ 

4'  ','  GUN  CALIBER  CALBER  =',G10.4,'  MM'/ 

5'  ','  BARREL  CROSS  SECTIONAL  AREA  A  =',G10.4,'  M**2'/ 

6'  ','  SPECIFIC  HEAT  RATIO  OF  PROPELLANT  GAS  GAM  =',G10.4/ 

7’  ','  MEAN  MOLECULAR  WEIGHT  OF  PROPELLANT  GAS  MBAR  =',G10.4, 

8  / 

9'  FLAME  TEMPERATURE  TF  =',G10.4,'  DEG  K'/ 

1'  ','  COVOLUME  ETA  =',G1 0.4,'  M**3/KG'/ 

2'  CHAMBER  VOLUME  CVO  =',G10.4,'  M**3'/) 

WRITE(6,6001) 

GO  TO  6000 

1  WRITE (IWRITE, 1000)  TITLE 

WRITE ( IWRITE, 1110)  V0 , L, MP , W, CALBER , A, GAM, MBAR , TA, ETA, 

1  CVO 

1110  FORMAT('0','  MUZZLE  VELOCITY  V0  =',G10.4,'  M/SEC'/ 

1'  ','  BORE  LEMiTO  L  r',G10.4,'  METERS'/ 

2'  ','  PROPELLANT  MASS  MP  =',G10.4,'  KG'/ 

3'  PROJECTILE  MASS  W  =',G10.4,'  KG'/ 

4'  ','  GUN  CALIBER  CALBER  =',G10.4,'  MM'/ 

5'  ','  BARREL  CROSS  SECTIONAL  AREA  A  =',G10.4,'  M**2'/ 

6'  SPECIFIC  HEAT  RATIO  OF  PROPELLANT  GAS  GAM  r',G10.4/ 

7'  ','  MEAN  MOLECULAR  WEIGHT  OF  PROPELLANT  GAS  MBAR  =',G10.4, 

8  / 

9'  AVERAGE  BARREL  GAS  TEMPERATURE  TA  =',G10.4,'  DEG  K»/ 

1'  ','  COVOLUME  ETA  =',G1 0.4,'  M**3/KG'/ 

2*  ','  CHAMBER  VOLUME  CVO  =' ,G10.4, '  M**3'/) 

WRITE (6, 6001) 

6000  CONTINUE 

CALCULATED  PARAMETERS 


L=L+CVO/A 
ALMP  =  A*L/MP 
LAMM=MP/W 
R=8317./MBAR 

TEFF  =  TF  -  . 5*(GAM~1.)»V0»V0*(. 3333+1 ./LAMM)/R 

IF(TA.NE.0.0)TEFF=TA 

EPS=ETA/(ALMP-ETA) 

M0=V0/SQRT(GAM»R*TEFF) 


INITIAL  CONDITIONS:  THET0 


C-  *1  -i-Eps 

G10=C+LAMM/(3.*GAM) 


36 


ooo  oooooo  oooo 


G1 0P=G1 0-LAMM* (C/M0+ . 5 )/GAM 
C 

J  =  0 

YOLD  =  (2.-M0+C)/(M0+C) 

EXP  =  2./(GAM+l.) 

20  J  =  J  +  1 

IF(J.GE.15)  GO  TO  20 

DENOM  =  M0  +  1.  +  EPS  +  (M0-1 .-EPS)/Y0LD 
YNEW  =  (EPS  +  (2.*G10/DENOM)**EXP)  /  G10 
C 

IF(ABS(YNEW-YOLD).LE.. 00001)  GO  TO  30 
YOLD  =  YNEW 
GO  TO  20 
C 

30  THET0  =  YNEW 


INITIAL  CONDITIONS:  THET0P 

ALPHA=  .5*(GAMf  1 . )/( 1  .-EPS/('mET0*G10) ) 

BETA=-1 . *THET0* ( M0+C ) / ( M0-C) 

BPRIME=-1 .+LAMM/(GAM*M0*M0) 

CPRIME-  -GAM')-GAM*FPS 

DELTA=(BPRIME*(1 .+1HET0)+CPRIME*(THET0-1 . )/M0)/( 1 .-C/M0) 
■IHET0P  =THET0«((1  .+THET0)*(1  .-G10P/G10)+DELTA+ 

1 G1 0P*( 1 . -ALPHA) * (BETA-1 . )/G1 0 )/( 2 .+ALPHA* (BETA-1 . ) ) 


SUBROUTINE  ZIHETA  INTEGRATES  DIFFERENTIAL  EQ.  FOR  IHETA 
ZIHETA  CALLS  RUNGE-KUTTA  SCHEME  SUBROUTINES  GILL  AND  GILL1 
ZIHETA  CALLS  DERIV  TO  CALCULATE  DERIVATIVES  OF  Z  AND  IHETA 

CALL  ZIHETA( DELTAU , DLTAU0 ,  IHET0 ,  IHET0P , IMAX, M, I , G1 0 , EJ ECTO) 

PRINT  MUZZLE  PROPERTIES 


DO  60  K=1,I 

T=TTAU(K)  »  L/V0  *  1000. 

RH02rEPS/(ETA«H1(K)) 

V2=SQRT(GAM*R*TEFF)*H1  (K)/(H1  (K)-EPS)**(  .5*(GAMi-1 . ) ) 
T2=TEFF*(H1 (K)-EPS)**(1 .-GAM) 

P2rR»T2/( 1 ./RH02-ETA) 

P2ATM  =  P2  /  101325. 

SONICV  =  SQRT(  GAM*R«T2  ) 

MACHNO  =  V2  /  SONICV 

FEJECT(K)  =FEJECT(K)/  M0  /  (1 .-ETA*MP/(A*L)) 
IF(MACHNO.LT.1.0)GO  TO  1001 
ZSQ= .476*GAM*P2ATM 
ZSQ1=SQRT(ZSQ) 

XM2=(  .4*ZSQ)**(GAM-1 .  )*(GAM4-1 . )/( (GAM-1 .  )**GAM) 

DO  2000  J=1 ,50 
AB=1 ./(GAM-1.) 

ABC=AB*GAM 

ABCD=(2.+(GAM-1.)*XM2) 

ABCDE=2. *GAM*XM2- (GAM-1 . ) 
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XFMP=-  .49*GAM**2»(  1  ./(GAMi-1 . )  )**AB*(ABCD**AB-2.*ABCD»*ABC/ABCDE)/ 
1ABCDE 

XFM=ZSQ-.49*GAM*( 1 ./(GAMf  1 . ) )**AB*ABCD**ABC/ABCDE 
C  WRITE(1,2001)XFM,XM2 

IF(ABS(XFM).LT..001)  GO  TO  2002 
XM2=XM2-XFM/XFMP 
IF(J.EQ.50)STOP  10 
2000  CONTINUE 
2002  XM1=SQRT(XM2) 

PM0=R*TEFF*( 1 .0-LAMM/3.0)/(ALMP-ETA) 

PM0=PM0/101325. 

TM0=TEFF*(1 .0-(GAM-1 .0)»LAMM/(GAM*3.0)) 

GAMP1=GAM+1 .0 
GAMM1=GAM-1 .0 
AXX=GAMM1/GAM 

MR= (GAMP1 *P2ATM*»AXX-2.0)/GAMM1 
MR=SQRT(MR) 

TR=T2*GAMP1/(2.0+GAMMl*MR*»2) 

UR=SQRT(GAM*R*TR)»MR 

TNI =GAMP1/ ( 2 .0+GAMM1 »XM1 »*2) 

TN2=1 .0+2.0*GAMM1»(GAM»XM2+1 .0)*(XM2-1 .0)/(XM2*GAMPl»*2) 
TNrTN1*TN2«T2 

MNP2= ( 2 .0+GAMM1 *XM2 )/ ( 2 .0*GAM*XM2-GAMM1 ) 

UN=SQRT(MNP2»R*TN*GAM) 

EPSD=0 .69*SQRT(GAM*P2ATM) 

XXY=GAMP1/(2.0*GAMM1 ) 

XXYY= (GAMP1/ ( 2.0+GAMM1 *XM2 ) ) **XXY 
TETTA=0 . 96*EPSD**2«XM1 *XXYY 
UB=(  1 .0-TETTA)*UR+TET;rA*UN 
TB1=(1 .0~TETTA)*TETTA*GAMMl/2.0 
1B2=MR»*2*(1 .0-(UN/UR)»*2)*TR 
TB=( 1 .0-TETTA)*TR+TETTA»TN+TBl*TB2 
PR=1 .0 

UMS=SQRT(GAM*R»T2) 

AB=P2ATM*(TB/T2)*(UMS/UB)»A 

RB=SQRT(AB/PI) 

WRITE(6,5000) 

5000  F0RMAT(32X,48HMUZZLE  GAS  PROPERTIES  WHEN  PROJECTILE  IS  EJECTED,/) 
WRITE (6, 6001) 

WRITE(6,5001)PM0 

5001  FORMAT(40X,9HPRESSURE=,2X,F7.2,4H  ATM) 

WRITE(6,5002)TM0 

5002  FORMAT(40X,12HTEMPERA1URE=,2X,F7.1,2H  K) 

WRITE(6,5003)V0 

5003  FORMAT(40X,9HVELOCITY=,2X,F7.1,6H  M/SEC) 

WRITE(6,6001) 

6001  F0RMAT(2X,///) 

WRITE(6,5004) 

5004  F0RMAT(32X,56HMUZZLE  GAS  PROPERTIES  WHEN  MUZZLE  VELOCITY  BECOMES  S 
IONIC,/) 

WRITE(6,6001) 

WRITE(6,5001)P2ATM 

WRITE(6,5002)T2 

WRITE(6,5003)UMS 

FEJ=FEJECT(K) 
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WRITE(6,5005)FEJ 

5005  FORMAT(il0X,31HFRACTION  OF  EJECTED  PROPELLANT=,2X,F6.4) 

WRITE (6, 6001) 

WRITE(6,5006) 

5006  FORMAT(32X,69HFLCW  CONDITIONS  AT  REFLECTED  SHOCK  WHEN  MUZZLE  VELOC 
IITY  BECOMES  SONIC,/) 

WRITE (6, 6001) 

WRITE(6,5001)PR 

WRITE(6,5002)TR 

WRITE(6,5003)UR 

ALPHU1.0-TETTA 

WRITE(6,5007)ALPH1 

5007  FORMAT(40X,41HFRACTION  OF  GAS  ENTERING  REFLECTED  SHOCKr,2X,F6.4) 
WRITE(6,6001) 

WRITE(6,5008) 

5008  FORMAT (32X,66HFLCW  CONDITIONS  AT  NORMAL  SHOCK  WHEN  MUZZLE  VELOCIH 
1  BECOMES  SONIC,/) 

WRITE(6,6001) 

WRITE(6,5001)PR 

WRITE(6,5002)TN 

WRITE(6,5003)UN 

WRITE(6,6001) 

WRITE (6, 5009) 

5009  FORMAT(32X,76HFLOW  CONDITIONS  AT  MIXING  REGION  BOUNDARY  WHEN  MUZZL 
IE  VELOCITY  BECOMES  SONIC,/) 

WRITE(6,6001) 

WRITE(6,5001)PR 

WRITE(6,5002)TB 

WRITE(6,5003)UB 

WRITE(6,5010)RB 

5010  FORMAT(40X,17HBOUNDARY  RADIUS  =,2X,F6.3,'  M' ,////) 

C 

WRITE(  9 ,501 1  )TN ,  TM0 ,  TB,  PM0 ,  UB,  ALPI11 ,  RB,  XMAX,  PRNT,  FDL,  KEY 

5011  FORMAT(10F8.3,I1) 

C  THUS  IS  DATA  PASSED  ON  TO  MTOB,  BLAKE,  CONCEN,  AND  LAPP 
C 

5555  FORMATCIH  ,''MACH  NO  IS  .GE.  1") 

IFCMACHNO.GE.I .0)WRITE(6,5555) 

IF(MACHNO.GE.1.0)STOP 

C 

1001  CONTINUE 
C 

WRITE( IWRITE, 1090)  T, RH02 , V2 , P2 , P2ATM, T2 ,SONICV, MACHNO, FEJECT(K) 
1090  FORMATS  ',^013. 4) 

C 

60  CONTINUE 

IFd.GT.DSTOP 


SUBROUTINE  RATEAU  CALCULATES  THE  MUZZLE  PROPERTIES  AFTER 
THE  RAREFRACTION  WAVE  HITS  THE  BREECH 


WRITE(IWRITE,1100) 

1100  FORMAT('0  RAREFRACTION  WAVE  HAS  HH  THE  BREECH*//) 
T  =  TTAU(I)*L/V0*1000. 
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T2  =  TEFF*(H1(I)-EPS)**(1 .-GAM) 

EJECTO  =  EJECTO/M0/(1.-ETA*MP/(A»L)) 

CALL  RATEAU(T,T2, EJECTO, I, JJ) 

C 

C  PLANE  A 
C 

WRITE(IWRITE,2110) 

2110  FORMATCI ',T45,' PLANE  A  FLOW  PROPERTIES'//*  TIME',T10, 

1  'DENSITY', T20, 'TEMP', T30, 'VELOCITY', T43,'XI/D',T50, 'SHOCK', 

2  T60, 'PLANE  A' ,T70, 'C02* ,T80, 'CO' ,T90, *H20' ,T100, 'H2'/ 

3*  MSEC,  T11,'KG/M3',T21,*  DBG  K*  ,T31 , 'M/SEC  ,T50,*  LENGTH' , 

4  T60 , '  DIAMETER ' ,  T70 ,  ’  M0L/M3 ' ,  T80 ,  *  M0L/M3 ' ,  T90 ,  *  M0L/M3 ' , 

5  T100,'MOL/M3'/) 

C 

DO  70  1=1, JJ 

P0  =  PE(I)»((GAM+1.)*.5)**(GAM/(GAM-1.)) 

XID  =  .69*SQRT(GAM»PE(I)) 

XID  =  XID/(1.  +  .197*PHI**.65) 

SHOCKL  =  1 .25*XID 

HAREA  =  A*SQRT(.5*(GAM-1.))*(2./(GAMf1.))»*((GAM+1.)*.5/ 

1  (GAM-1.)) 

HAREA  =  HAREA/SQRT(1.-P0**((1.-GAM)/GAM)) 

HAREA  =  HAREA*P0**(1./GAM) 

HDIA  =  SQRT(HAREA/A) 

RHOA  =  RHOE(I)/PE(I)**(1./GAM) 

TA  =  TE(I)/PE(I)**((GAM-1 .)/GAM) 

UA  =  UE(I)»(A/HAREA)»PE(I)**(1./GAM) 

WRITE (IWRITE, 1120)  TIME( I ) , RHOA, TA, UA, XID , SHOCKL, HDIA, 

1  C02(I),C0(I),H20(I),H2(I) 

1120  FORMAT(1P11G10.3) 

70  CONTINUE 

ISENTROPIC  REGION  AND  MACH  DISC 

WRITE( IWRITE, 1130) 

1130  FORMATCI  *,T45,' ISENTROPIC  REGION  FLOW  FIELD'// 

1 15 ,  'TIME' ,  120, '  R/D' , T35 , ' MACH  NO. ' , T50, '  RHO'  ,T65, 

2  'PRESSURE' ,T80, 'TEMPERATURE' ,T95, ' VELOOCITY*/ 

3T6 , ' MSEC* , T50 , ' KG/M3 ' , T66 , *  ATM' , T8l , ' DBG  K* , T96 , ' M/SEC ' /) 

C 

DO  90  KK=5,JJ,5 

XID  =  .69*SQRT(GAM«PE(KK)) 

XID  =  XID/(1.  +  .197*PHI**.65) 

STEP  =  .1*XID 

ROVERD  =  .69*SQRT(2.*GAM)  -  STEP 
80  ROVERD  =  ROVERD  +  STEP 

IF(ROVERD.GE.XID)  ROVERD  =  XID 
CALL  MACHN(ROVERD,FMACH,PHI,GAM,XID) 

RHOI  =  RH0E(KK)*((1.+GAM)/(2.+(GAM-1.)»FMACH»FMACH))»* 

1(1. /(GAM-1.)) 

RATIO  =  RHOI/RHOE(KK) 

PI  =  PE(KK)*( RATIO )**G AM 
TI  =  TE(KK)*(RATIO)**(GAM-1.) 

UI  =  UE(KK)»FMACH»(RATIO)**(.5*(GAM-1.)) 

WRITEdWRITE,  1140)  TIME (KK) , ROVERD,  FMACH, RHOI,  PI, TI, UI 
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1140  F0RMAT(1P7G15.4) 

IF(ROVERD.LT.XID)  GO  TO  80 
FMACHS  =  FMACH*FMACH 

FMACH2  =  (1.  +  .5*(GAM-1.)*FMACHS)/(,GAM*FMACHS-.5*(GAM-1.)) 
FMACH2  r  SQRT(FMACH2) 

RH0MD2  =  RH0I*(GAMf1.)*FMACHS/(( GAM-1. )*FMACHS+2.) 

PMD2  =  1. 

TMD2  =  TE(KK)*(2.*GAM*FMACHS-GAMf1.)/(GAMf1.)/FMACHS 
UMD2  =  RH0I/RH0MD2*UI 

WRITE ( IWRITE, 1150)  TIME ( KK) , ROVERD,  FMACH2 , RH0MD2 , PMD2 , TMD2 , UMD2 
1150  FORMAT(1P7G15.4//) 

90  CONTINUE 
C 

WRITE(IWRITE,1160) 

1160  FORMAT('0»,'  THIS  IS  IHE  END  OF  IHE  PROGRAM* »S  OUTPUT*) 

C 

STOP 

END 


SUBROUTINE  ZTHETA(DELTAU,DLTAU0 ,1HET0,'IHET0P, IMAX,M, I, 

1  G10,EJECTO) 

ZIHETA  DETERMINES  Z(TAU)  AND  IHETA(TAU)  USING 
RUNGE-KUTTA  INTEGRATION  SCHEME 

IMAXrDIMENSION  OF  Z  IN  MAIN  PROGRAM 
I=HIGHEST  SUBSCRIPT  FOR  WHICH  VALUES  STORED  IN  Z 
M=.NUMBER  OF  ITERATIONS  BETWEEN  STORRED  VALUES 
DOUBLE  Y(20),DERY(20),Q(20),G1,G1P 
REAL  M0,LAMM 
EXTERNAL  DERIV 
COMMON  /G/  Q,A(4),B(4),C(4) 

COMMON/WORKB/THETA(500) ,Z(500) , HI (500) ,H1P(500) ,IHETAP(500) , 
1  TTAU(500)  ,FEJECT(500) 

COMMON/WORKA/LAMM, GAM, EPS , M0 , B0 , C0 , D0 , KPMAX 

INITIAL  CONDITIONS 

THETAC 1 ) =THET0+THET0P*DLTAU0 
TAU=DLTAU0 
TTAU(l)  =  DLTAU0 
Y(1)=THETA(1) 

DERY(l)  =  THET0P 
EXPON  =  -.5*(GAM+1.) 

DERIV  DETERMINES  DY/DX  AND  THE  VALUES  OF  G1  AND  G1P 

CALL  DDERIV 
NSUB=1 

CALL  DERIV(TAU,Y,DERY,G1 ,G1P,NSUB) 

C  CALL  DERIV(TAU,Y,DERY,G1,G1P) 

NSUB=2 

CALL  DERIV(TAU,Y,DERY,G1,G1P,NSUB) 

Z(1)=  Y(2) 
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IHETAP(l)  =  DERY(I) 

H1(1)=Y(1)*G1 

H1P(1)=DERY(1)»G1+Y(1)*G1P 

FEJECT(l)  =  .5*(((THET0*G10-EPS)**EXPON)  +((Y(1)*G1-EPS) 
1**EXPON))*DLTAU0 
EJECT  =  FEJECTd) 

G10LD  =  G1 


GILL  INITIALIZES  RUNGE  -  KUTTA  ROUTINE  GILL1 


CALL  GILL(1,Q,DERIV) 

CALL  GILL( 1 ) 

1=1 

10  1=1+1 

IF(I.GE.IMAX)  GO  TO  50 

GILL1  INTEGRATES  ONE  STEP  REPLACING  TAU  AND  Y  WITH  NEW  VALUES 

Y10LD  =  THETA(I-I) 

Y20LD=Z(I-1) 


DO  30  J=1,M 

CALL  GILL1 (DELTAU , TAU , Y, DERY) 

STORING  CONDITION:  Z(TAU0)=0  ;  LINEAR  INTERPOLATION  DETERMINES 
TAU0  &  THETA (TAU0) 

IF(Y(2).GE.0.0)  GO  TO  20 
R=Y20LDA(2) 

TAU0=TAU-DELTAU/( 1 .-R) 

TTAU(I)  =  TAU0 
Y(2)  =  0. 

Z(I)  =  Y(2) 

C 

SLOPE  =  (Y(1)-Y10LD)/DELTAU 
Y(1)  =  SLOPE* (TAU0-TAU)  +  Y(1) 

THETA(I)=Y(1) 

C 

C  CALL  DERIV(TAU0,Y, DERY, G1,G1P) 

NSUB=2 

CALL  DERIV ( TAU0 , Y,  DERY,  G1 ,  G1 P ,  NSUB ) 

THETAP(I)  =  DERY(I) 

H1(I)=Y(1)*G1 

H1P(I)=DERY(1)*G1+Y(1)*G1P 

FEJECT(I)  =  EJECT+.5*(((Y10LD*G1QLD-EPS)**EXP0N)+ 

1 ( ( Y ( 1 ) *G1 -EPS ) **EXPON ) ) ♦DELTAU 
EJECTO  =  .5*((Y(1)*G1-EPS)**EXPON)*DELTAU 
RETURN 

C  20  CALL  DERIV  ( TAU ,  Y,  DERY,  G1 ,  G1 P ) 

20  NSUB=2 

CALL  DERIV ( TAU, Y, DERY, G1 ,G1P,NSUB) 

EJECT  =  EJECT  +.5*(((Y10LD*G10LD-EPS)**EXP0N)+ 

1 ( ( Y ( 1 ) *G1 -EPS ) ♦*EXPON ) ) «DELTAU 
G10LD  =  G1 
Y10LD  =  Y(1) 
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Y20LD  =  Y(2) 

30  CONTINUE 
C 

raETA(I)=Y(l) 

TTAU(I)  =  TAU 
Z(I)  =Y(2) 

IHETAP(I)  =  DERY(I) 

H1(I)=Y(1)*G1 

H1P(I)=DERY(1)*G1+Y(1)*G1P 
FEJECT-(I)  r  EJECT 
C 

GO  TO  10 
C 

50  RETURN 
END 

SUBROUTINE  DERIV(TAU,Y,DERY,G1 ,G1P,NSUB) 
SUBROUTINE  DDERIV 

DERIV  CALCULATES  DY/DX  Y(1)=THETA  DERY(1 )=IHETP 

Y(2)=Z  DERY(2)=ZP 

DOUBLE  Y(2) ,DERY(2) ,G1 ,G1P,B,  C,F,F1 ,F2 
REAL  M0,LAMM 

COMMON  /WORKA/LAMM, GAM, EPS , M0 , B0 , C0 , D0 , KPMAX 
IF(NSUB.NE.I)  GO  TO  30 
CONSTANTS 

AA=LAMM/(GAM*M0*M0) 

A1=AA/(2.-GAM) 

A2=l .+EPS-LAMM/(6.»GAM)+AA/(GAM-1 .) 

A3=A1/(1.-GAM) 

A4=LAMM/(2.«GAM) 

A8rl .+EPS-AA*.5 
A9=-1.*EPS*(GAMi.1.)*.5 
A10=(1.-GAM)*A1*.5 

BB  =  GAM  -  1 . 

BBB  =  2.  -  GAM 

B1  =  1.  +  AA/BB  -  2./(BB»M0)*(l.+AA/(3.*BB)) 

B2  =  2./(BB*M0) 

B3  =  1.  -  AA/2. 

B4  =  -.5»AA*BB/BBB 
B5  =  -.5*EPS»BB 

B6  =  (GAM  +  1.  )»AA/(6.*BB*BBB) 

B7  =  -1.*AA/BBB 
B8  =  AA/(BB*BBB) 

RETURN 

CALL  DERIV(TAU,Y,DERY,G1,G1P) 

30  CONTINUE 
C  EQUATIONS 
C 

F=1 .+TAU 
C 

FI  =  B3  +  (B4+B5mu)/F  +  B6/F«*BB 
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F2  =  1 .  +  B7/F  +  B8/F**BB 
Y(2)  =  F  *  (B1  +  B2/F**(.5*BB)  *  F1)/F2 
C 

B=1 ./F-A1/(F*F)+A1/F»*GAM 

C=F**( .5*( 1 .-GAM) )*(A8+A9*TAU/F+A10/F+A1 *.5/F**(GAM-1 . ) ) 
DERY(2)=B*Y(2)-C/M0 
C 

G1=A1+A2»F+A3*F**(2.-GAM)+A4*Y(2)*Y(2)/F 
G1P=A2+(2.-GAM)*A3*F**(1 .-GAM)+A4*(-1.*Y(2)«Y(2)/(F*F)+ 
12.*Y(2)*DERY(2)/F) 

IF(TAU.EQ.0.)  GO  TO  60 
FACT=Y(1)*G1-EPS 
C  WRITE(6,*)FACT 

IF( FACT . LT. 0 . ) FACT= 1 . E-50 

DERY(1)=-1.*Y(1)*(1.+Y(1))«G1P/G1+Y(1)/(1.-Y(2))»(C*(1.-Y(1))/M0 
1 -B*Y( 2)*( 1 .+Y( 1 ) ) +2. *Y( 1 ) *01 *FACT**(- . 5*(GAM+1 . ) ) /M0) 

C 

C 

60  CONTINUE 
RETURN 
END 

SUBROUTINE  GILL(N) 

SUBROUTINE  GILL(N,Q,DERIV) 

DOUBLE  Y(N),DERY(N),Q(N),G1,G1P,XK,DY 
DOUBLE  Y(20),DERY(20),Q(20),G1,G1P,XK,DY 
COMMON  /G/  Q,A(4),B(4),C(4) 

DATA  A(1),B(1),C(1),C(4)/0.5,-1.0,-0.5,-0.5/ 

DO  10  1=1, N 
10  Q(I)=0.0 

C1=1.0/SQRT(2.0) 

A(2)=1.0-C1 
A(3)=1.0+C1 
A(4)=1 .0/6.0 
B(2)=-A(2) 

B(3)=-A(3) 

B(4)=-1.0/3.0 
C(2)=B(2) 

C(3)=B(3) 

RETURN 
END 

SUBROUTINE  GILL1 (DX, X, Y, DERY) 

C  CALL  GILL1(DX,X,Y, DERY) 

DOUBLE  Y(20) ,DERY(20) ,Q(20) ,G1 ,G1P,XK,DY 
COMMON  /G/  Q,A(4),B(4),C(4) 

DO  30  1=1,4 

IF(I.BQ.2.0R.I.EQ.4)  X=X+0.5*DX 
C  CALL  DERIV  (X, Y,DERY,G1 ,G1P) 

NSUB=2 

CALL  DERIV(X,Y,DERY,G1,G1P,NSUB) 

DO  30  J=1,N 
XK=DX*DERY(J) 

DY=A(I)*XK+B(I)*Q(J) 

Y(J)=Y(J)+DY 

Q(J)=Q(J)+3.0*DY+C(I)*XK 
30  CONTINUE 
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RETURN 

END 

SUBROUTINE  RATEAU(TSTAR , T20 , EJECTO, II, J) 

REAL  M0,LAMM,L,MP,MACHNO 
COMMON/WORKA/LAMM, GAM, EPS , M0 , B0 , C0 , D0 , KPMAX 
COMMON/WORKB/THETAC  500 ) , Z ( 500 ) , H1( 500 ) , HI P( 500) , THETAP( 500 ) , 
1  TTAU(500),FEJECT(500) 

COMMONAJORKC/TEFF, ETA, A, L, MP , V0, ALMP , R 
COMMON/WORKD/PE(50) ,UE(50) ,TE(50) ,RHOE(50) ,TIME(50) , 

1  TEMP(50),  FKP(50), 002(50), CO(50),H2O(50),H2(50) 

DATA  DELTAT/1.E-03/ 

DATA  IWRITE/6/ 

TSTAR  =  TSTAR/1000. 

IHET  =  (.5»(GAMf1.))*»((GAMf1.)/(GAM-l.))/GAM 

THET  =  'mET/R/TEFF/(l.+LAMM/6.)/(1.-ETA*LAMM)**(GAM-1.) 

IHET  =  2.*L*(1.+.013*ETA/ALMP)/(GAM-1.)»SQRT(THET) 

T0  =  ((1.+LAMM/6.)*TEFF/T20)**(1. /(GAM-1.)) 

10  =  T0*(ALMP-ETA)/(ALMP+1.07*ETA) 

T0  =  T0  +  2.07*ETA/(ALMP+1.07*ETA) 

10  =  THET*(T0**(.5*(GAM-1.))  -  1.) 

T  =  TSTAR  -  DELTAT 
EJECT  =  FEJECT(II) 

J  =  0 

10  J  =  J  +  1 

T  =  T  +  DELTAT 

TMSEC  =  1000. »T 

TREPLC  =  T  -  TSTAR  +  T0 

RH02  =  (1.  +  TREPLC/THET)**(2./(GAM-1.)) 

RH02  =  (ALMP  +  1 .07*ETA)*RHO2  -  1 .07*ETA 

RH02  =  1./RH02 

P2  =  R*TEFF*(1.+LAMM/6.) 

P2  =  P2*(ALMP-ETA)**(GAM-1.) 

P2  =  P2/(1 ./RH02-ETA)**GAM 
T2  =  P2»(1./RH02-ETA)/R 
V2  =  (ALMP-ETA)/(1./RH02-ETA) 

V2  =  V0/M0/(1.-ETA*RHO2)*V2**(.5*(GAM-1.)) 

SONICV  =  SQRT(GAM*R«T2) 

MACHNO  =  V2/S0NICV 
P2ATM  =  P2/101330. 

PE(J)  =  P2ATM 
UE(J)  =  V2 
TE(J)  =  12 
RHOE(J)  =  RH02 
TIME(J)  =  TMSEC 

DETERMINE  THE  REACTION  RATE  CONSTANT  FKP 

TA  =  T2/P2ATM*»( (GAM-1. )/GAM) 

DO  20  1=1,  KPMAX 
IF(TA.GE.TEMP(I))  GO  TO  20 
IF(I.LE.I)  GO  TO  30 

SLOPE  =  (FKP(I)  -  FKP(I-1))/(TEMP(I)  -  TEMP(I-I)) 

FKPP  =  SLOPE*(TA  -  TEMP(I))  +  FKP(I) 

GO  TO  40 
20  CONTINUE 
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30  WRITE(IWRITE,1000)  TMSEC,TA 

1000  FORMAT(‘  EXIT  PLANE  TEMPERATURE  T2  IS  OUTSIDE  THE  RATE', 

1  '  CONSTANT  INTERPOLATION  TABLE  RANGE'/ •  COMPUTATION', 

2  '  TERMINATED  Wm  TIME  =',F10.3,'  ,T2  =',F10.3/) 

40  CONTINUE 

AA  =  (D0*FKPP  +  B0  +  C0)**2  +  4.»B0*C0*(FKPP-1 . ) 

AA  =  SQRT(AA)  -  D0*FKPP  -  B0  -  C0 
AA  =  AA/(2.*(FKPP-1.)) 

C02(J)  =  AA 
CO(J)  =  B0  -  AA 
H20(J)  =  C0  -  AA 
H2(J)  =  D0  +  AA 

IFCJ.EQ.D  GO  TO  50 
EJECTN  =  .5WMP*RH02»V2»DELTAT 
EJECT  =  EJECT  +  EJECTN  +  EJECTO 
EJECTO  EJECTN 

50  WRITE(IWRITE,1010)  TMSEC,RH02,V2,P2,P2ATM,T2,S0NICV,MACHN0,EJECT 
1010  FORMATC  ',1P9G13.4) 

IF(P2ATM.LE.1.)  RETURN 

GO  TO  10 

END 

SUBROUTINE  MACHN(R,FM,PHI,GAM,XI) 

MACHN  CALCULATES  THE  MACH  NUMBER  GIVEN  THE  AXIAL  DIATANCE 
FROM  THE  MUZZLE  BY  HALF  INTERVAL  SEARCH 

DATA  IWRITE/6/ 

DATA  EPS/.01/ 

F(X)  =  SQRT(.49*GAM*(2.+(GAM-1.)»X*X)**(GAM/(GAM-1.))/ 

1  (1.  +  .197*PHI**.65)**2/(2.*GAM*X»X  -  GAM  +  1.)/ 

2  (GAM^1.)*♦(1./( GAM-1.))) 

FMIN  =  1. 

FMAX  =  1. 

DO  10  1=1,20 
FMAX  =  FMAX  +  FMAX*FMAX 
IF(F(FMAX).GE.XI)  GO  TO  30 
10  CONTINUE 
DO  20  1=1,200 
FM  =  .5*(FMIN  +  FMAX) 

DELTA  =  F(FM)  -  R 
IF(ABS(DELTA).LE.EPS)  RETURN 
IF(DELTA.LE.0.)  FMIN  =  FM 
IF(DELTA.GT.0.)  FMAX  =  FM 
20  CONTINUE 
30  CONTINUE 

WRITE(IWRITE,1000)  FMIN, FMAX,  FM,R, DELTA 
1000  FORMAT( '0WARNING-SUBROUTINE  MACHN"S  CALCULATION', 

1  '  OF  MACH  NUMBER  HAS  NOT  CONVERGED  AFTER  200  ITERATIONS'/ 

2'  FMIN  =',1PG10.3,'  FMAX  =',G10.3,'  FM  =',G10.3,'  R  =', 

3  G10.3,'  DELTA  =',G10.3) 

STOP 

END 
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ooo  o  oo  oo  oo  ooooo 


PROGRAM  MrOB( INPUT, OUTPUT, TAPE9 , TAPE4 , TAPES , TAPE6=0UTPUT) 

TAPES  IS  OUTPUT  TAPE  FOR  BLAKE,  CONCEN,  AND  LAPP 
TAPE9  IS  INPUT  TAPE  FROM  MEFF 

TAPE4  IS  INPUT  BOILERPLATE  FOR  IHE  APPROPRIATE  GAS 

DIMENSION  A(30,S0),TITLE(20) 

REWIND  4 
REWIND  S 
REWIND  9 

READ  IHE  RESULTS  OF  THE  MEFF  CALCULATION 
READ(9,400)TITLE 

READ( 9 , 300) T1 , T2 , T3 , P, V, ALPHA, RADIUS, XMAX, PRNT, FDL, KEY 

NOW  READ  THE  BLAKE  BOILERPLATE 
DO  20  1=1,30 
IF(EOF(4))30,20 
20  READ(4,100)(A(I,J),J=1,S0) 

HERE  ONCE  BOILERPLATE  FILE  IS  READ 
30  11=1-2 

DO  40  K=1,II 

40  WRITE(S,100)(A(K,J),J=1,S0) 

NOW  FOR  THE  SPECIAL  LINES 

50  WRlfE(S,120) 

120  FORMAT(9HPRL,CON,2) 

WRITE(S,110)T1 

110  FORMAT(11HPOI,P,1.,T,,FS.3) 

WRITE(S,130)P,T2 

130  FORMAT(6HPOI, P, ,FS.3,3H,T, ,FS.3) 

WRITE(S,170) 

170  FORMAT(4HQUIT) 

WRITE(S,1S0) ALPHA 
1S0  FORMAT(FS.3) 

WRITE(S,400)  TITLE 
WRITE(S , 1 90)RADIUS, XMAX, PRNT 
190  FORMAT(3E10.3) 

WRITE(S,210)FDL,KEY 
210  FORMAT(FS.3,I1) 

WRITE(S,200)T3,V 
200  FORMAT(2F8.3,/,4H*EOI) 

300  FORMAT(10FS.3,I1) 

100  FORMAT(S0A1) 

400  FORMAT(20A4) 

END 
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TIT,M30A1 

PRL,C0N,2 

REJ ,  N,  K2S04 ,  C,  C2,  CH,  CH20,  HN03 
REJ,C(S),K2S04$ 

REJ,K0H$,K02$,K202$ 

REJ, H2S, S20, S02 , K$, K20, K202 
REJ , KCO$ , KSO$ , K20$ , K$ 

REJ,K2C03$ 

REJ,K2S$ 

UNI, ENG 

CM2,NC1260,27.9,NG,22.42,NQ,46.84,EC,1 .49,KS,1 ., 
ALC,.25,C,.1 
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I 


PROGRAM  QDNCEN (INPUT, OUTPUT, TAPE1 ,TAPE2,TAPE5=INPUT, 

+  TAPE6r OUTPUT, TAPES) 

C 

C  13  SPECIES 
C 

DIMENSION  NAM(17) ,CONC(17,2) ,NA(2,29) ,CO(2,29) 

DIMENSION  CON(17),SUM(2) 

DATA  CONST/ 'CONST'/ 

C 

C  THE  NAME  OF  EACH  SPECIES  MUST  BE  4  CHARACTERS 
C 

DATA  NAM/'  H20' , '  CO','  H2','  N2','  C02','  H','  OH','  O', 

+'  02','  K','  KOH','  K02','  H02'/ 

REWIND  1 
REWIND  2 

READ (8, 300) ALPHA 
NSPEC=13 
K=0 
1=0 

SUM(1)=0. 

SUM(2)=0. 

1  READ (1,200) ALINE 
IF(ALINE.NE.CONST)GOTO  1 
1=1+1 

READ (1,200) JUNK 
DO  10  J=l,29 

READ(1,100)NA(I,J),CO(I,J) 

IF(NA(I,J).EQ.'  ')GOTO  20 
10  CONTINUE 
20  IFd.EQ.DGOTO  1 
50  K=K+1 

DO  30  I=1,NSPEC 
DO  40  J=1 ,29 

IF(NA(K,J).EQ.NAM(I))CONC(I,K)=CO(K,J) 

40  CONTINUE 
30  CONTINUE 

DO  60  L=1 ,NSPEC 
60  CONTINUE 

IF(K.EQ.1)G0T0  50 
DO  70  1=1,2 
DO  80  J=1,NSPEC 
SUM(I)=SUM(I)+CO(I,J) 

80  CONTINUE 
70  CONTINUE 
WRITE (6, 400) 

WRITE(6,500)ALPHA 
DO  90  I=1,NSPEC 

C0N(I)=(1.-ALPHA)*(C0NC(I,1)/SUM(1))+ALPHA*(C0NC(I,2)/SUM(2)) 

WRITE(6,400)NAM(I),CON(I) 

90  CONTINUE 

WRITE(2,700)(CON(I),I=1,NSPEC) 

WRITE (6, 600) 

STOP 

100  FORMAT(4X,A4,14X,1PE11.5) 

200  FORMAT(9X,A5) 
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300  FORMAT(F8.3) 

400  FORMAT(2X,A4,5X,1PE10.3) 

500  F0RMAT(1H1,2X,  •  ALPHA=  SFS. 
600  FORMAT(IHI) 

700  FORMAT(7E10.4) 

END 
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87.07 

47.535 

4000. 

13.905 

87.72 

50.313 

4200. 

13.925 

88.33 

53.097 

H02 

33. 

008 

0.5 

100. 

7.949 

61.574 

-1.596 

200. 

8.003 

55.132 

-0.799 

400. 

8.907 

54.717 

0.877 

600. 

9.980 

56.116 

2.771 

800. 

10.769 

57.657 

4.850 

1000. 

11.365 

59.123 

7.066 

1200. 

11.831 

60.481 

9.387 

1400. 

12.197 

61 .734 

11.791 

1600. 

12.485 

62.891 

14.261 

1800. 

12.714 

63.965 

16.782 

2000. 

12.895 

64.966 

19.343 

2200. 

13.041 

65.902 

21.937 

2400. 

13.160 

66.781 

24.558 

2600. 

13.256 

67.610 

27.200 

2800. 

13.336 

68.393 

29.859 

3000. 

13.403 

69.135 

32.534 

3200. 

13.459 

69.840 

35.220 

3400. 

13.507 

70.511 

37.917 

3600. 

13.547 

71.153 

40.622 

3800. 

13.582 

71 .766 

43.335 

4000. 

13.612 

72.354 

46.055 

4200. 

13.638 

72.918 

48.780 

CO 

+0 

+M  = 

C02  +M 

26 

7.0E-33  0. 

-4369. 

CO 

+02 

+  = 

C02  +0 

15 

4.2E-12  0. 

-47664. 
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0 

+0 

+M 

=02 

+M 

25  3.0E-34 

0. 

1792.0 

CO 

+0H 

+ 

=C02 

+H 

16  2.8E-17- 

■1.3 

660.0 

OH 

+H2 

+ 

=H20 

+H 

161.90E-15- 

-1.3 

-3626. 

H 

+02 

+ 

rOH 

+0 

152.40E-10 

0. 

-16393. 

0 

+H2 

=0H 

+H 

16  3.0E-14 

-1. 

-8902. 

OH 

+0H 

=H20 

+0 

151.05E-11 

0. 

-1093. 

H 

+H 

+M 

=H2 

+M 

223.00E-30 

1. 

0. 

H 

+0H 

4.M 

=H20 

+M 

231.00E-25 

2. 

0. 

H 

+02 

+M 

=H02 

+M 

25  1.5E-32 

0. 

994. 

H 

+H02 

+ 

=0H 

+0H 

15  1.7E-10 

0. 

-994. 

CO 

+H02 

+ 

=C02 

+0H 

15  2.5E-10 

0. 

-23645. 

H 

+H02 

+ 

=H2 

+02 

15  4.2E-11 

0. 

-695. 

H 

+H02 

+ 

=H20 

+0 

15  8.5E-12 

0. 

-994. 

OH 

+H02 

+ 

=H20 

+02 

11  3.0E-11 

0. 

0. 

0 

+H02 

+ 

=0H 

+02 

11  3.5E-11 

0. 

0. 

0 

+H 

+M 

=0H 

+M 

22  1.0E-29 

1. 

0. 

H02 

+H2 

+ 

=H20 

+0H 

15  1.0E-12 

0. 

-18678. 

H 

+KOH 

+ 

=H20 

+K 

15  1.8E-11 

0. 

-1987. 

K 

+0H 

+M 

=KOH 

+M 

22  1.5E-27 

1. 

0. 

K02 

+0H 

+ 

=KOH 

+02 

11  2.0E-11 

0. 

0. 

K 

+02 

+M 

=K02 

+M 

22  3.0E-30 

1. 

0. 

K 

+H02 

+ 

=K02 

+H 

15  1.0E-11 

0. 

-13000. 

K02 

+H2 

+ 

=KOH 

+0H 

15  3.0E-12 

0. 

-19870. 

APPENDIX  J 
LAPP  PROGRAM  LISTING 
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c 

C  RELEASE  VERSION  OF  LAPP  MDDIFIED  TO  ACCEPT  UP  IHROUGH  49  REACTIONS. 

C  MARCH  1983 

C 

C 

PROGRAM  LAPP ( INPUT, OUTPUT, TAPES , TAPES , TAPE5=INPUT, TAPE6=OUTPUT, 

+  TAPE?, TAPES, TAPE10,TAPE11,TAPE1 2) 

C 

C  TAPE2  IS  INPUT  FROM  PROGRAM  CONCEN 
C  TAPES  IS  INPUT  FROM  FILE  OF  LAPP  DATA 
C  TAPES  IS  INPUT  FROM  PROGRAM  MTOB 
C 

c*»*»*»*»** 

c**********  AXISTMMETRIC  MIXING  WITH  NON-EQUILIBRIUM  CHEMISTRY  ** 

C**********  AEROCHEM  RESEARCH  LABORATORIES  ** 

C**********  PRINCETON,  N.  J.  ** 

C«**»***»** 

c********************^****************************^*******************^* 

DIMENSION  A(30),RHO(30),Y(30),T(30),PSI(30),RT(30),SUM(30),AR(25), 
1HSTAT(30),H(25,30),ALPHA(25,30),RALPHA(25,30),CP(25,30),SIGMA(30), 
2  WTMOLE(25),CPBAR(30),C(25,9),AID(25),ETA(30),RATIO(30), 

3RU(30),U(30),TITLE(12),XLE(30),XMU(30)  ,G(25) ,WTMIX(30) 

4RC(49,3),IRRR(49,5),WP(25),WM(25),WDOT(25,30),SAVET(30),SAVEU(30), 

5  IRR(49),FREQ(30),SAVEA(25,30),  PC(4) ,ZID(5) , 

6  ECC(30) ,HOUT(30) , YOUT(30) , RHOOUT(30) ,XMUOUT(30) , XLT(30) , 
7T4(30),TFDG(30),IRT(49),RP(49,30),RM(49,30) 

DIMENSION  XAME(6),  ISAVE(6),  FRBQA(6),  ALOC(50,6) 

DIMENSION  1^(30,24)  ,CPTB(25,30)  ,HTB(25,30)  ,GTB(25,30)  ,HF(25) 

DIMENSION  TABLE(16),ZSPEC(16) 

COMMON/CONSTS/AMULT, AMULT1 , AMULT2 , AMULT3 
COMMON/UNITS/NUNITA, NUNITB, NUNITC, NUNITD, NUNITE, NUNITF, 

*NUNITG, NUNIIH , NUNITI, NUNITJ , NUNITK, NUNITL , NUNITM, NUNIIN , 
*NOUT,NDBG,NNNOUT 

COMMON/A/  CM(25,25,26),CM1(25,25),QX(25,26),QX1(25) 

COMMON/C/  IZSPEC,ISPEC(16) 


COMMON 

A 

f 

RHO 

7 

Y 

) 

T 

7 

PSI 

7 

RT 

COMMON 

SUM 

1 

AR 

7 

HSTAT 

7 

H 

7 

ALPHA 

7 

RALPHA 

COMI'ION 

CP 

f 

SIGMA 

7 

WIMOLE 

7 

CPBAR 

7 

C 

COMMON 

AID 

J 

ETA 

7 

RATIO 

7 

RU 

7 

U 

1 

TITLE 

COMMON 

XLE 

7 

XMU 

7 

G 

7 

WTMIX 

7 

WDOT 

COMMON 

SAVEU 

7 

SAVET 

7 

WM 

f 

WP 

7 

RC 

COMMON 

SAVEA 

7 

PC 

7 

X 

f 

XMAX 

COMMON 

PRNT 

7 

DXMIN 

7 

DX 

7 

FDL 

7 

DELPSI 

7 

RJ 

COMMON 

XK2 

7 

P 

7 

ZID 

7 

FREQ 

7 

ECC 

7 

DPDX 

COMMON 

YOUT 

7 

HOUT 

7 

RHOOUT 

7 

IRRR 

7 

IRR 

• 

IFINIS 

COMMON 

IPAGE 

7 

MPSI 

7 

MY 

7 

NS 

7 

NR 

f 

7 

lEDGE 

COMMON 

IIURB 

7 

IPRESS 

7 

NPSI 

f 

ITEST 

7 

ITER 

• 

lECC 

COMMON 

IRT 

7 

XMUOUT 

7 

XLT 

1 

T4 

7 

TFDG 

1 

lOUT 

COMMON 

IOUT1 

7 

I0UT2 

,RP 

7 

RM 

7 

ISAVE 

7 

IPUNCH 

COMM3N  TKINET, NFREQA, ALOC, FREQA, QQ1 00 , QQ200 , QQ300 , QQ400 

c 

C»***  SETTING  XAME  FOR  IHOSE  SPECIES  INVOLVED  IN  COLLISION  FREQ.  CALC. 
C 

DATA  XAME(l)/6HCO  / 


67 


DATA  XAME(2)/6HC02 

.  / 

DATA  XAME(3)/6HH20 

/ 

DATA  XAME(4)/6HH2 

/ 

DATA  XAME(5)/6HN2 

/ 

DATA  XAME(6)/6HHCL 

/ 

DATA  TABLE/8HCO 

,8HHCL 

,8HH20 

,8HOH 

8HC02 

,8HHF 

,8HCN 

,8HO 

'  8HBF 

,8HBFO 

,8HBH02 

,8HBF2 

8HBF3 

,8HBO 

,8H 

,8HAL203 

CCCCC  COMMENT  OUT  NAMELIST  DEFINITION  WHICH  FOLLOWS  IF  PROCESS  NAMELIST 
CCCCC  BY  CALLS  TO  SUBROUTINE  NMLST. . . .ALSO  COMMENT  OUT  NORMAL  READS 
CCCCC  AND  WRITES  FOR  IBM  SUPPORTED  NAMELIST 

CCCCC  NAMELIST/NUNITS/NUNITA, NUNITB, NUNITC, NUNITD, NUNITE, NUNITF, NUNITG, 
CCCCC*NUNI1H , NUNITI , NUNITJ , NUNITK, NUNITL , NUNITM , NOUT, NDBG 
C 

C  INITIALIZE  VARIABLES 
AMULT=0.3M8 
AMULT3=14.5939 
JPR0C=2 
NUNITA=2 
NUNITB=3 
NUNITC=5 
NUNITD=7 
NUNITE=7 
NUNITF=8 
NUNITG=0 
NUNITH=0 
NUNITI=0 
NUNITJ=0 
NUNITK=0 
NUNITL=0 
NUNITM=0 
NUNITN=5 
NOUT=6 
IINOUT=6 
NNNOUT=10 
NDBG=6 

C 

AV=6.025E23 

C...  SET  CONVERSION  TO  ENGLISH  UNITS  MULTIPLIERS 
AMULT1=AMULT«AMULT 
AMULT2=AMULT*AMULT1 
AMULT3=SQRT(AMULT3) 

4  IFINIS=0 
DO  9999  1=1,8536 
9999  A(I)=0.0 

DO  300  I  =  1,6 
300  ISAVE(I)  r  0 
R=82.06 

C .  BEGIN  TO  INPUT  FILE  DESIGNATED  BY  NUNITB 

READ  (NUNITF,333)(TI1LE(I),I=1,10) 

IFCNUNITM. LT.0 )WRITE(NNOUT, 333)  (TITLE( I ) , 1=1 , 1 0) 

C 

C  NS=  NO.  OF  SPECIES 
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C  NR=  NO.  OF  REACTIONS 
C 

c 

READ( NUNITB, 1 00)  MPSI, NS, ITURB, NR, lOUTl , IOUI2 , IPUNCH, ITIME, 
»IPRESS,NT 

IF(NUNITM.LT.0)WRITE(NNOUT,100)  MPSI,NS,I'njRB,NR,IOUT1 ,10UT2, 
•IPUNCH, ITIME, IPRESS,  NT 
READ  (NUNITB, 111)  (FRBQA(I),  1=1,6) 

IF(NUNITM.LT.0)  WRITE (NNOUT,l 11)  (FREQA(I),I=1 ,6) 

DO  113  1=1,6 

IF  (FREQA(I))  113,  114,  113 

113  CONTINUE 
1=7 

114  NFREQA  =  I-l 
NPSI=MPSI-1 

READ  ( NUNITB, 1 000 ) X, XMAX, PRNT, XLE ( 1 ) , SIGMA( 1 ) , RJ , XK2 
READCNUNITF, 1 000)RJ, XMAX, PRNT 

IF(NUNITM. LT.0 )  WRITE(NNOUT, 1000)X, XMAX, PRNT, XLE( 1 ) , SIGMA( 1 ) , RJ , XK 
X=X/AMULT 
XMAXrXMAX/AMULT 
PRNT=PRNT/AMULT 
RJ=RJ/AMULT 
DX=0.1*RJ 
C 

C  INPUT  MINIMUM  STEPSIZE  LIMIT  (DXMIN) 

C 

READ(NUNITB,555)  DXMIN, FDL,PC(1),PC(2),PC(3),PC(4),lHOT,TCOOL 
READ( NUNITF, 556 )FDL, KEY 
556  FORMAT(F8.3,I1) 

IF(NUNITM.LT.0)  WRITE(NNOUT,555)DXMIN,FDL,  (PC(KI)  ,KI=1 ,4)  ,‘IHOT, 
•TCOOL 

IFCKEY.EQ.DGOTO  98765 

NOUT=10 

NNOUT=10 

NCBG=10 

NNNOUT=6 

98765  DXMIN=DXMIN/AMULT 
PC(2)=PC(2)*AMULT 
PC(3)=PC(3)*AMULT1 
PC(4)=PC(4)*AMULT2 

READ  (NUNITB, 1 000) P, T( 1 ) , T(MPSI ) , U( 1 ) , U(MPSI ) , DELPSI, TKINET 
READ( NUNITF, 987)T( 1 ) ,U( 1 ) 

987  FORMAT(2F8.3) 

IF(NUNITM.LT.0)  WRITE(NNOUT, 1000 )P,T(1 ) ,T( MPSI) ,U( 1 ) ,U(MPSI) , 
•DELPSI, TKINET 
U(1)=U(1)/AMULT 
U(MPSI)=U(MPSI)/AMULT 
DELPSI=DELPSI/AMULT3 
IF(TKINET.EQ.0.0)  TKINET  =  400.0 
C 

C***«  the  value  of  30  SECONDS  IS  TO  ALLOW  FOR  COMPILE  TIME 
C 

ILIMIT  =  60^ITIME 
IDIFFT  =  0 
CALL  TICK(ISECST) 
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IF((ISECST+60*ITIME).GT. 86400)  IDIFFT  =  86400-ISECST 
UNIT  =  U(1) 

IQ77  =  2 
USUB01  =0.0 
C 

C  TURBULENCE  MODELS 
C 

IF  (ITURB  -  3)  8600,9010,9010 
8600  IF  (ITURB  -  1)  9011,9010,9011 

9010  USUB01  =  0.95*  (U(l)-U(MPSI))  +  U(MPSI) 

9011  CONTINUE 

IF(DELPSI)  3011,3012,3011 
C 

C»»««  READING  OF  CENTERLINE  CONCENTRATIONS  FROM  FILE  PRODUCED 
by  BLAKE  AND  CONCEN. 

C 

3012  READ  (NUNITA,1000)(ALPHA(J,1),J=1,NS) 

C 

C  BE  SURE  NO  DENSITIES  ARE  ZERO 
DO  7777  UK=1,NS 

ALPHAdJK,  1  )=AMAX1  (ALPHACIJK,  1 ) ,  1  .E-99) 

7777  CONTINUE 
C 

IF(NUNITM.LT.0)  WRITE(NNOUT, 1000) (ALPHA( J, 1 ) , J=1 ,NS) 

C 

C»»«»  NOW  READ  THE  CONCENTRATIONS  ON  IHE  EDGES  -  AMBIENT  AIR. 
C 

READ(NUNITB,1000)(ALPHA(J,MPSI),J=1,NS) 

IF( NUNITM.lt. 0)  WRITE(NNOUT,1000)(ALPHA(J,MPSI),J=1,NS) 

MM0D=MPSI-2 

DO  4001  I=1,MM0D 

T(I)=T(1) 

U(I)=U(1) 

DO  4001  J=1,NS 

4001  ALPHA(J,I)=ALPHA(J,1) 

DO  4002  J=1,NS 

4002  ALPHA(J,NPSI)=ALPHA(J,MPSI) 

GO  TO  3015 

3011  READ  (NUNITB, 1000) (T(I), 1=1, MPSI) 

IF(NUNITM.LT.0)  WRITE(NNOUT, 1000) (T(I) , 1=1 , MPSI) 

READ  (NUNITB, 1000) (U(I), 1=1, MPSI) 

IF( NUNITM. LT. 0 )  WRITE( NNOUT , 1 000) ( U ( I ) , 1=1 , MPSI ) 

DO  7  1=1, MPSI 
U(I)=U(I)/AMULT 

READ  (NUNITB,1000)(ALPHA(J,I),J=1,NS) 

IF(NUNITM.LT.0)  WRITE(NNOUT,1000)(ALPHA( J,I), J=1 ,NS) 

7  CONTINUE 
C 

C  NEW  THERMO  DATA  DATA  INPUT  IN  JANNAF  TABLE  FORM 
C 

3015  DO  1991  1=1, NS 

READ( NUNITB, 222)  AID(I) ,WTMOLE(I) ,HF(I) 

IF( NUNITM. LT.0)  WRITE(NN0UT,222)  AID(I) ,WTMOLE(I),HF(I) 

DO  10  IT=1,  NT, 2 

ITP1=IT+1 
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READ(NUNITB,102)  TTB(IT,I),  CPTB(I, IT) ,GTB(I,IT) IT) , 
»TTB(ITP1 ,1) ,CPTB(I, ITP1 )  ,GTB(I, ITPl ) ,HTB(I, ITP1 ) 

IF(NUNITM.LT.0)  WRITE(NNOUT,102)TTB(IT, I) ,CPTB(I, IT) jGIBd, IT) , 
*HTB(I, IT) , TTB( ITPl , I) , CPTB( I, ITPl ) ,GTB( I, ITPl ) , HTB( I, ITPl ) 
G1B(I,IT)=  -GIB(I,IT)*TTB(IT,I)  +HF(I)  *1000. 

GTB(I,ITP1)=  -GTB(I,ITP1)*TTB(ITP1,I)  +HF(I)*1000. 
HTB(I,IT)=(HIB(I,IT)+HF(I))*1000. 

HTBd,  ITPl  )=(HTBd,  ITPl )+  HFd)  )*1000. 

10  CONTINUE 

IF(WTMOLE(I)-1.0)  1972,1991,1991 
1972  lECCrl 
1991  CONTINUE 

MODIFICATIONS  FOR  LAPP/ARC  INTERFACE  PROGRAM 

WRITE(NUNITD,333)  (TITLEd),I=l  ,10) 

IDENTIFY  LAPP  SPECIES  IN  IHE  ARCTABLE 

TABLErSPECIES  NAMES  AS  FOUND  IN  ARCTABLE 
15=NUMBER  OF  SPECIES  IN  ARCTABLE 
IZSPEC=NUMBER  OF  LAPP  SPECIES  FOUND  IN  ARCTABLE 


M=0 

DO  3  1=1, NS 
DO  3  J=1,16 

IF(AIDd).NE.TABLE(J))  GO  TO  3 
M=M+1 

ZSPEC(M)=AIDd) 

ISPEC(M)=I 
3  CONTINUE 
IZSPECrM 

WRITE(NUNITD,334)  IZSPEC 
334  FORMAT(8I10) 

WRITE (NUNITD, 333)  (ZSPEC(M)  ,M=1  ,nSPEC) 

WRITE  (NUNITD, 102)  THOT,TCOOL 
WRITE(NUNITD,1000)  RJ 
DO  301  J  =  1 ,6 
DO  301  I  =  1,NS 

IF(AIDd).EQ.XAME(J))  ISAVE(J)  =  I 
301  CONTINUE 

DO  1992  1=1, NR 

READ  (NUNITB,444)  (ZID(J) ,  J=1 ,5)  ,IRRd)  ,IRTd) ,  (RCd,K)  ,K=1 ,3) 
IF(NUNITM.LT.0)  WRITE(NNOUT,444) (ZID( J) , J=1 ,5) ,  IRRd) ,  IRTCD , 
*(RCd,K),K=l,3) 

DO  1993  J=1,5 
IRRRd,J)=0 
DO  1993  L=1,NS 

IF(ZID(J)-AID(L))  1993,1994,1993 
1994  IRRRd,J)=L 
1993  CONTINUE 
1992  CONTINUE 

DO  912  I=1,MPSI 

W'IVR=0.0 

DO  632  J=1,NS 
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632  WTVR=WTVR+ALPHA(J,I)*WTMOLE(J) 

DO  633  J=1,NS 

633  ALPHA(J,I)=ALPHA(J,I)/WTVR 
912  CONTINUE 

IF(DELPSI)  903,3041,903 
3041  DUM=0.0 

DO  5001  J=1,NS 
5001  DUM=DUM+ALPHA(J,1) 

XMD=MMOD~1 

DELPSI=SQRT(P«U( 1 )/42.285E0/T( 1 )/DUM)*RJ/XMD 
903  DO  20  1=1 ,29 
XI=I-1 

PSI(I)=XI»DELPSI 

XLE(I)=XLE(1) 

20  SIGMA(I)=SIGMA(1) 

DO  90  I=NPSI,29 
RT(I)=T(MPSI) 

T(I)=T(MPSI) 

DO  80  J=1,NS 

RALPHACJ, I)=ALPHA( J,MPSI) 

80  ALPHA(J,I)=ALPHA(J,MPSI) 

RU(I)=U(MPSI) 

90  U(I)=U(MPSI) 

IF(NOUT.GT.0)CALL  INOUT 
PPUNCH  =  P 
P=2117.0*P 
DPDX=0.0 
C 

C  PRESSURE  OPTION,  IF  IPRESSr  0,  PRESSURE 
C  IS  CONSTANT  AND  =  TO  P,  IF  IPRESS  =  1 ,  COEFFICIENTS  CALLED 

C  PC(1),PC(2),PC(3),PC(4)  ARE  INPUT. 

C  EQUATION  USED  IS  P  =  PC(1)  +  PC(2)*X  +  PC(3)*X*X  +  PC(4)*X*X*X 

C 

2  IF(IPRESS)  821,822,821 

821  P=(PC( 1 )+X*(PC(2)+X»(PC(3)+X*PC(4) ) ) )*21 17.0 
DPDX=(PC(2)+X»(2.0*PC(3)+X*3.0*PC(4)))»2117.0 
PPUNCH  =  P/2117.0 

822  DO  31  I=1,MPSI 
WTMIX(I)=0.0 
DO  30  J=1,NS 

30  WTMIX(I)=WTMIX(I)+ALPHA(J,I) 

RHO(I)=P/89517.501/T(I)/WTMIX(I) 

31  RHOOUT(I)=RHO(I)/1.94 
DO  805  I=1,MPSI 

C 

C  FREE  STREAM  VELOCITY  WILL  BE  SET  TO  1 .0  FPS  IF  ZERO  IS  ENTERED 
C 

U(I)  =  AMAX1(1.0E0,U(I)) 

TFDG( I ) =AMAX1 ( T( I ) , 1 00 . 00E0 ) 

T4(I)=TFDG(I)**4 

XLT(I)=ALOG(TFDG(I)) 

CPBAR(I)=0.0 

HSTAT(I)=0.0 

TX=T(I) 

DO  805  J=1,NS 
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JJ=J 

CALL  TKEY  (TX,TTB,ITKEY,SDT,Hiyr,NT,  JJ) 

IF  (ITKEY.EQ.0)  GO  TO  9  »  »  ' 

CP(J,I)=0.0 

H(J,I)=0.0 

CALL  LIPLN( ITKEY, J, CPTB, SDT, HDT, AX) 

CP(J,I)=AX»45055.31 

CALL  LIPLNC ITKEY, J,  HTB,SDT,HDT,AX) 

H(J,I)=AX*45055.31 

HSTAT(I)=HSTAT(I)+H(J,I)»ALPHA(J,I) 

CPB^(I^=CPBAR(I)+CP(J,I)»ALPHA(J,I) 

Y(1)=0r0* 

C********CALL 

ETA(2)=DELPSI/SQRT(RH0( 1 )*U(2) ) 

Y ( 2) =DELPSI/SQRT ( RHO ( 2 ) *U( 2 ) ) 

DO  25  I=3,MPSI 

1003  FORMAT(I5,7E10.3) 

IF(TEMP,LT.0.0)  CALL  OUTPUT 
Y(I)  =SQRT(TEMP) 

25  CONTINUE 
C 

C****  HAS  MIXING  REGION  INTERSECTED  X  AXIS  YET,  YES  IF  0  OR  - 

IF  (ITURB-e)  8010,8011,8010 

8010  IF  (U(1)  -  USUB01)  9000,9000,9001 

C  MODEL  6  COMMON  CALCULATIONS 

c 

8011  QQ1  =  (U(1)+U(MPSI))/2.0 
DO  8012  I=2,MPSI 

8012  8013,8013,8012 

8013  QQ2  =  (QQ1-U(I-1))/(U(I)-U(I-1)) 

QQ100  =  Y(I-1)+(Y(I)-Y(I-1))»QQ2 
QQ30  =  T(I-1)+(T(I)-T(I-1))*QQ2 

.  QQ3  =  CPBAR(I-1)+(CPBAR(I)-CPBAR(I-1))»QQ2 
^  =  1*0/(VJTMIX(I-1)+(WTMIX(I)-WTMIX(I-1))*QQ2) 

QQ5  =  89517. 501/QQ4 
QQ6  =  QQ3/(QQ3-QQ5) 

QQ7  =SQRT(QQ6*QQ5»QQ30) 

QQ300  =  QQ1/QQ7 
QQ8  =  (ABS(U(1)-U(MPSI)))/2.0 
_  IF  (QQ300-1.2)  8014,8014,8015 

8014  ^00  =^.0468+QQ300*((QQ300»(-.0460))+.0256*(QQ300*QQ300)))*XK2 

8015  QQ400  =  .0248*XK2 

^16  IF  ((U(1)-USUB01)»(U(1)-U(MPSI)))  8020  8020  8800 

8020  IF  (USUB01 -9000,0)  8021,8021,8900  ’  ’ 

8021  USUB01  =  10000.0  ’ 
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XDUM=X»0.30M8 
WRITE(NOUT,9900)  XDUM 
GO  TO  8900 
C 

C  MODEL  6  BEFORE  MIXING  ZONE  REACHES  AXIS 
C 

8800  QQ9  =  0.95*  (U(I)-U(MPSI))  +  U(MPSI) 

DO  8802  I  =  2,MPSI 

IF  ((QQ9-U(I))*(QQ9-U(I-1)))  8804,8804,8802 

8802  CONTINUE  ^  x  x  ^  ^ 

8804  QQ200  =  Y(I-1)+(Y(I)-Y(I-1))*(QQ9-U(I-1))/(U(I)-U(I-1)) 

QQ10  =  QQ400*QQ8*(QQ100-QQ200) 

DO  8810  I  =  1,MPSI 
8810  XMU(I)  =  QQ10*RHO(I) 

GO  TO  98 
C 

C  MODEL  6  AFTER  MIXING  ZONE  REACHES  AXIS 
^8900  QQ11  =  QQ400*QQ100*QQ8 

DO  8910  I  =  1,MPSI 
8910  XMU(I)  =  QQ11*RHO(I) 

QQ200  =0.0 
GO  TO  98 

9000  1077  =  9 
USUB01  =0.0 
WRITE(NOUT,9900)X 

9001  LL  =  ITURB  +  IQ77 
C 

C  EDDY  VISCOSITY  MODELS 

^  GOTO  (91 ,99, 8666, 78, 8667, 8668, 9003, 91 ,99, 45, 78, 26, 33, 78), LL 

C 

C  MODEL  1  BEFORE  MinNG  ZONE  REACHES  AXIS 

^8666  XMU( 1 ) =0 .001 37*(X+1 .0E-05) *ABS(RHO( 1) *U ( 1 ) -RHO(MPSI ) *U (MPSI ) ) 
GO  TO  37 
C 

C  MODEL  3  BEFORE  MIXING  ZONE  REACHES  AXIS 

8667  XMU ( 1 ) =0 . 00 1 37 * ( X+1 . 0E-05 ) *RHO ( 1 ) *ABS( U ( 1 ) -U ( MPS I ) ) 

GO  TO  37 
C 

C  MODEL  4  BEFDRE  MIXING  ZONE  REACHES  AXIS 

^866 8  XMU ( 1 ) =0 . 00 1 37 * ( X+1 . 0E-05 ) »RHO (MPSI ) *ABS( U( 1 ) -U ( MPSI ) ) 

GO  TO  37 

91  DO  92  1=1 ,MPSI 
C 

C  MODEL  0  LAMINAR  FLOW 

c 

92  XMU(I)=3.05E-8*T(I)*«1.5/(T(I)+111.0) 

GO  TO  98 

45  DUM= .5»(RHO( 1)*U( 1 )+RHO(MPSI)*U(MPSI)) 

DO  52  J=1,MPSI 
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I=MPSI-J+1 

IF(RHO(I)*U(I)-DUM)  52,52,51 
52  CONTINUE 

MODEL  1  AFTER  MIXING  ZONE  REACHES  AXIS 
51  Z=Y(I)+(Y(I)-Y(I+1))*(RH0(I)*U(I)-DUM)/(RH0(I)»U(I)-RH0(I+1)*U(I+1 

D) 

XMU ( 1 ) =XK2*Z*ABS( RHO( 1 ) *U ( 1 ) -RHO(MPSI ) *U (MPSI ) ) *0.025 
GO  TO  37 

99  DO  39  I=1,MPSI 
39  XMU(I)=XK2*0.025 
GO  TO  98 

78  RD=(U(1)+U(MPSI))/2.0 
DO  47  I=2,MPSI 

IF  ((RD-U(I))*(RD-U(I-1)))  48,48,47 

47  CONTINUE 

48  RHALVE=ETA(I-1)+(ETA(I)-ETA(I-1))*(RD-U(I-1))/(U(I)-U(I-1)) 
DUMMY=XK2*RHALVE*ABS( U( 1 )-U(MPSI) ) *0.025 
XMU(1)=DUMMY*RH0(1) 

DO  79  I=2,MPSI 

MODEL  2  BEFORE  MIXING  ZONE  REACHES  AHS 

79  XMU(I)=DUMMY*(RH0(1)*ETA(I)/Y(I))**2/RH0(I) 

GO  TO  98 

MODEL  2  AFTER  MIXING  ZONE  REACHES  AXIS 

26  RD=(U(1)+U(MPSI))/2.0 
DO  27  I=2,MPSI 

IF  •  ( (RD-U( I) ) *(RD-U(I-1 ) ) )  28 ,28,27 

27  CONTINUE 

28  RHALVE=Y(I-1 )+(Y(I)-Y(I-1 ) )*(RD-U(I-1 ) )/(U(I)-U(I-1 ) ) 

MODEL  3  AFTER  MIXING  ZONE  REACHES  AXIS 

XMU ( 1 ) =XK2*RHALVE*RH0( 1 ) *ABS( U( 1 ) -U(MPSI) ) *0 .025 
DO  29  I=1,MPSI 

29  XMU(I)=XMU(1) 

GO  TO  98 

33  RD=(U(1)+U(MPSI))/2.0 
DO  34  Ir2,MPSI 

IF  ((RD-U(I))*(RD-U(I-1)))  35,35,34 

34  CONTINUE 

35  RHALVE=Y(I-1 )+(Y(I)-Y(I-1 ) )*(RD-U(I-1 ) )/(U(I)-U(I-1 ) ) 

MODEL  4  AFTER  MinNG  ZONE  REACHES  AXIS 

XMU( 1 )=XK2*RHALVE*RH0(MPSI)*ABS(U( 1 )~U(MPSI) )*0.025 
37  DO  36  I=1,MPSI 

36  XMU(I)=XMU(1) 

GO  TO  98 

MODEL  5  BEFORE  MIHNG  ZONE  REACHES  AXIS 
9003  XMU(l)  =  0.00137*(X+1.0E-05)*ABS(UNIT-U(MPSI))*RHO(1) 
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DO  9004  I  =  2,MPSI 

XMU(I)=0.00137*(X+1.0E-05)*ABS(UNIT-U(MPSI))*(RHO(1)**2/RHO(I)) 
9004  CONTINUE 
98  A(1)=0.0 
C 

C  CALCULATE  A 
C 

DO  44  I=2,MPSI 

44  A(I)r:XMU(I)*RHO(I)*U(I)*Y(I)»Y(I)/PSI(I) 

DO  899  L=1,NPSI 
RRT=1.986«fT(L) 

ROOTT=SQRT(T(L)) 

TX=T(L) 

DO  855  1=1, NS 
11=  I 

CALL  TKEYdX, TIB,  ITKEY,  SDT,  HDT,  NT,  II) 

IF  (ITKEY. EQ.0)  GO  TO  9 
G(I)=0.0 
WP(I)=0.0 
WM(I)=0.0 
QX(I,L)=0.0 
DO  872  J=1,NS 
872  CM(I,J,L)=0.0 

CALL  LIPLN(ITKEY,I,  GIB,  SDT, HDT, AX) 

G(I)=AX 
855  CONTINUE 
C 

C  REACTION  CALCULATION 

C  REACTION  KINETICS  CONTINUE  DOWN  TO  400  DEGREES  K 
C  UNLESS  TKINET  IS  SET  TO  A  VALUE  OTHER  THAN  400  K 

C  REACTION  KINETICS  FOR  ALL  REACTIONS  CONTINUE  DOWN  TO  TKINET 
C 

IF(T(L)-.TKINET)  3256,3256,3259 
3259  CONTINUE 

DO  862  1=1, NR 
RP(I,L)=0.0 
RM(I,L)=0.0 
KK  =  IRT(I) 

C 

C  REACTION  CONSTANT  TYPE 
C 

GO  TO  (841, 842, 843, 844, 845, 846, 847), KK 

841  RATE=RC(I,1)«AV 
GO  TO  849 

842  RATE=RC(I,1)/T(L)«AV 
GO  TO  849 

843  RATE=RC(I,1)/T(L)/T(L)*AV 
GO  TO  849 

844  RATE=RC(I,l)/ROOTT*AV 
GO  TO  849 

845  RATE=RC(I,1)«EXP(RC(I,3)/RRT)*AV 
GO  TO  849 

846  RATE=RC(I,1)«EXP(RC(I,3)/RRT)/T(L)*»RC(I,2)*AV 
GO  TO  849 

847  RATE=RC(I,1)/T(L)/R00TT«AV 
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849  CONTINUE 
K=IRR(I) 

C 

C  TYPE  OF  REACTION 
C 

GO  10(864,865,866,870,871, 834 , 835 , 836 , 837 , 838 ) , K 

870  J1=IRRR(I,1) 

J2=IRRR(I,2) 

J3=IRRR(I,3) 

E  =  (G(J1)+G(J2)-G(J3))/RRT 
IF(ABS(E).LT.80.0)  GO  TO  700 
IF(E.LT.0.0)  E=EXP(-80.0E0) 

IF(E,GT.0.0)  E  =EXP(80.0E0) 

GO  TO  701 

700  E  =EXP(E) 

701  CONTINUE 

CRR=RATE*RHOOUT(L) 

RP(I,L)=CRR*RH00UT(L)*ALPHA(J1 ,L)*ALPHA(J2,L) 
RM(I,L)=CRR*ALPHA(J3,L)/E/R/T(L) 

DO  771  J=1 ,  3 
SIGN=1 .0 

IF  (J.GT.2)  SIGN=-1.0 
IROW=  IRRR(I,J) 

CM(IR0W,J1,L)=  CM(IRCW,J1,L)  +SIGN*RP(I,L)/ALPHA(J1 ,L) 
CM(IR0W,J2,L)=  CM(IR0W,J2,L)  +SIGN»RP(I,L)/ALPHA(J2,L) 
CM(IRCW,J3,L)=  CM(IR0W,J3,L)  -SIGN*RM(I,L)/ALPHA(J3,L) 

771  QX(IROW,L)=  QX(IROW,L)+  SIGN*RP(I,L) 

GO  TO  868 

871  J1=IRRR(I,1) 

J2=25 

J3=IRRR(I,3) 

J4=IRRR(I,4) 

E  =  (G(J1)-G(J3)-G(J4))/RRT 
IF(ABS(E).LT.80.0)  GO  TO  702 
IF(  E.LT.0.0)  E  =EXP(-80.0E0) 

IF(  E.GT.0.0)  E  =EXP(80.0E0) 

GO  TO  703 

702  E  =EXP(E) 

703  CONTINUE 

CRR=  RATE«RHOOUT( L) *RHOOUT( L)  »WTMIX( L) 
RP(I,L)=CRR*ALPHA(J1,L) 

RM( I , L) =CRR*R«T(L) »RHOOUT( L) *ALPHA( J3 , L) *ALPHA( J4 , L)/E 
DO  772  J=1 ,  4 
SIGN=1 .0 

IF  (J.GT.2)  SIGN=-1.0 
IROWr  IRRR(I,J) 

IF(J.BQ.2)  IR0W=25 

CMdROW,  J1  ,L)=  CMdRCW,  J1  ,L)+SIGN»CRR 

CM(  IROW,  J3 ,  L)=  CM(  IROW,  J3 ,  L)-SIGN*RMd,  L)/ALPHA(  J3 ,  L) 

CM ( IROW , J4 , L) =  CM ( IROW , J4 , L) -SIGN*RM( I , L) /ALPH A ( J4 , L) 

772  QX(IROW,L)=  QX(IROW,L)-  SIGN*RM(I,L) 

GO  TO  867 

864  J1=IRRR(I,1) 

J2=IRRR(I,2) 

J3=IRRR(I,3) 
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J4=IRRR(I,4) 

E  =  (G(J1)+G(J2)-G(J3)-G(J4))/RRT 
IF(ABS(E).LT.80.0)  go  to  704 
IF(E.LT.0.0)  E  =EXP(-80.0E0) 

IF(E.GT.0.0)  E  =EXP(80.0E0) 

GO  TO  705 

704  E  =  EXP(E) 

705  CONTINUE 

CRR=RATE*RHOOUT( L) »RHOOUT(L) 
RP(I,L)=CRR*ALPHA(J1,L)*ALPHA(J2,L) 

RM( I , L) =CRR*ALPHA( J3 , L) *ALPHA( J4 , L)/E 
DO  773  J=1 ,  4 
SIGN=1 .0 

IF  (J.GT.2)  SIGN=-1.0 
IROW=  IRRR(I,J) 

CM(IROW,J1,L)=  CM(IRCM,J1,L)  +SIGN*RP(I,L)/ALPHA( J1 ,L) 
CM(IR0W,J2,L)=  CM(IROW,J2,L)  +SIGN*RP(I,L)/ALPHA( J2,L) 
CM(IR0W,J3,L)=  CM(IROW,J3,L)  -SIGN*RM(I,L)/ALPHA( J3,L) 
CM( IROW, J4, L) =  CMC IROW, J4 , L)-SIGN*RM( I, L) /ALPHAC  J4 , L) 

773  QX(IROW,L)=  QX(IRCW,L)+  SIGN*(RP(I,L)-RM(I,L)) 

GO  TO  867 

865  J1=IRRR(I,1) 

J2=IRRR(I,2) 

J3=IRRR(I,3) 

E  =  (G(J1)+G(J2)-G(J3))/RRT 
IF(ABS(E).LT.80.0)  GO  TO  706 
IF(E.LT.0.0)  E=EXP(-80.0E0) 

IF(E.GT.0.0)  E=EXP(80.0E0) 

GO  TO  707 

706  E  =EXP(E) 

707  CONTINUE 

CRR= RATE*RHOOUT( L) *RHOOUT( L) *WTMIX( L) *AV 
RPC I, L)=CRR*RHOOUTC  L) *ALPHAC  J1 , L) *ALPHAC  J2 , L) 
RMCI,L)=CRR*ALPHACJ3,L)/CE»R*TCL)) 

DO  774  J=1,3 
SIGN=1 .0 

IF  C J.GT.2)  SIGN=-1.0 
IROW=  IRRRCI,J) 

CMCIROW,J1,L)=  CMCIRCW,J1,L)  +SIGN*RP C I, L) /ALPHAC J1 ,L) 
CMCIR0W,J2,L)=  CMCIROW,J2,L)  +SIGN*RPCI,L)/ALPHAC J2,L) 
CMCIROW,J3,L)=  CMCIRCW,J3,L)  -SIGN*RMCI,L)/ALPHAC J3,L) 

774  QXCIROW,L)=  QXCIROW,L)+  SIGN*CRPCI,L  )) 

GO  TO  868 

866  J1=IRRRCI,1) 

J2=IRRRCI,2) 

J3=IRRRCI,3) 

J4=IRRRCI,4) 

J5=IRRRCI,5) 

E=  CGCJ1)+GCJ2)-GCJ3)-GCJ4)-GCJ5))/RRT 
IFCABSCE).LT.80.0E0)  GO  TO  708 
IFCE.LT.0.0)  E=EXPC-80.0E0) 

IFCE.GT.0.0)  E  =EXPC80.0E0) 

GO  TO  709 

708  E  rEXPCE) 

709  CONTINUE 
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CRR=RATE*RHOOUT(L) »RHOOUT(L) 

RPC  I, L) =CRR*ALPHA( J1 , L) *ALPHA( J2 , L) 

RM( I , L) =CRR*ALPHA( J3 , L) »ALPHA( J4 , L) *ALPH A( J5 , L) *RHOOUT( L) »R*T( L) /E 
DO  775  J=1 ,  5 
SIGN=1 .0 

IF  (J.GT.2)  SIGN=-1.0 
IROW=  IRRR(I,J) 

CM(IROW,Jl,L)=  CM(IRCW,J1,L)  +SIGN*RP(I,L)/ALPHA(J1 ,L) 
CM(IR0W,J2,L)=  CM(IR0W,J2,L)  +SIGN»RP(I,L)/ALPHA( J2,L) 
CM(IROW,J3,L)=  CM(IRCW,J3,L)  -SIGN*RM(I,L)/ALPHA(J3,L) 

CMC IROW, J4 , L) =  CMC IROW, J4 , L) -SIGN*RMC I, L) /ALPHAC  J4 , L) 

CMC  IROW ,  J5 ,  L)  =  CMC  IROW,  J5 ,  L)  -SIGN^RMCI ,  D/ALPHAC  J5 ,  L) 

775  QXCIROW,L)=  QXCIR0W,L)+  SIGN*CRPCI,L)-2.»RMCI,L)) 

GO  TO  861 

837  J1=IRRRCI,1) 

J2=IRRRCI,2) 

J3=IRRRCI,3) 

CRR=RATE«RHOOUTCL) 

RPC I , L) =CRR*RHOOUTC  L) *ALPHAC  J1 , L) »ALPHAC  J2 , L) 

RMCI,L)=0.0 
DO  776  J=1 ,  3 
SIGN=1 ,0 

IF  CJ.GT,2)  SIGN=-1.0 
IROWr  IRRRCI,J) 

CMCIRCW,J1,L)=  CM(IROW,J1,L)  +SIGN*RPCI,L)/ALPHAC J1 ,L) 

CMC IROW, J2,L)=  CMCIROW,J2,L)  +SIGN*RP(I,L)/ALPHAC J2,L) 

776  QXCIRCW,L)=  QXCIRCIW,L)+  SIGN*  RPCI,L) 

GO  TO  868 

838  J1=IRRRCI,1) 

J2=25 

J3=IRRRCI,3)  . 

J4=IRRRCI,4) 

CRR=RATE*RHOOU  TC  L) »RHOOUT( L)  *WTMIXC  L) 

RPCI,L)=CRR*ALPHA(J1,L) 

RMCI,L)=0.0 
DO  777  J=1 ,  4 
SIGN=1 .0 

IF  CJ.GT.2)  SIGN=-1.0 
IROW=  IRRRCI,J) 

IF  CJ.EQ.2)  IROW=25 

777  CMC IROW, J1,L)=  CMC IROW, J1 ,L)  +SIGN*CRR 
GO  TO  867 

834  J1=IRRRCI,1) 

J2=IRRRCI,2) 

J3=IRRRCI,3) 

J4=IRRRCI,4) 

CRR=RATE«RHOOUTCL) *RHOOUTCL) 

RP  C I, L) =CRR*ALPHAC  J1 , L) *ALPHAC  J2 , L) 

RMCI,L)=0.0 
DO  778  J=1 ,  4 
SIGN=1 .0 

IF  CJ.GT.2)  SIGN=-1.0 
IROWr  IRRRCI,J) 

CMC IROW, J1,L)=  CM(IRCW,J1,L)  +SIGN«RP(I,L)/ALPHACJ1 ,L) 

CMC IROW, J2,L)=  CMCIR0W,J2,L)  +SIGN*RPCI,L)/ALPHAC J2,L) 
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778  QXdRW^Dr  OX(IRWd)+  SIGN*  RP(I,L) 

GO  TO  867 

835  J1=IRRR(IJ) 

J2=IRRR(I,2) 

J3=IRRR(I,3) 

CRR=RATE*RHOOUT( L) «RHOOUT( L) *WTMIX( L) «AV 
RP (I, L) =CRR*RHOOUT( L) *ALPHA( J1 , L) *ALPHA( J2 , L) 

RM(I,L)=0-0 
DO  779  J=1 ,  3 
SIGNrI .0 

IF  (J.GT.2)  SIGN="1.0 
IROW=  IRRR(I,J) 

CM(IROW,J1,L)=  CM(IRCIW,J1,L)  +SIGN*RP(I,L)/ALPHA(J1  ,L) 
CM(IR0W,J2,L)=  CM(IR0W,J2,L)  +SIGN*RP(I,L)/ALPHA( J2,L) 

779  QX(IRCW,L)=  QX(IRCW,L)-«-  SIGN*  RP(I,L) 

GO  TO  868 

836  J1=IRRR(I,1) 

J2=IRRR(I,2) 

J3=IRRR(I,3) 

J4=IRRR(I,4) 

J5=IRRR(I,5) 

CRR=RATE*RHOOUT( L) *RHOOUT( L) 

RP ( I, L) =CRR*ALPHA( J1 , L) *ALPHA( J2, L) 

RM(I,L)=0.0 
DO  780  J=1 ,  5 
SIGNrI .0 

IF  (J.GT.2)  SIGN==.1.0 
IROW=  IRRR(I,J) 

CM(IRCW,J1 ,L)=  CM(IRCM,J1,L)  +SIGN*RP(I,L)/ALPHA(J1 ,L) 
CM(IROW,J2,L)=  CM(IR0W,J2,L)  +SIGN*RP(I,L)/ALPHA( J2,L) 

780  QX(IROW,L)=  QX(IRCM,L)+  SIGN*  RP(I,L) 

CALCULATE  WDOT 

861  WP(J5)=WP(J5)+RP(I,L) 

WM(J5)=WM(J5)+RM(I,L) 

867  WP(J4)=WP(J4)-i-RP(I,L) 

WM(J4)=WM(J4)+RM(I,L) 

868  WP(J3)=WP(J3)+RP(I,L) 

WM(J3)“WM(J3)+RM(I,L) 

WP(J2)=WP(J2)+RM(IjL) 

WM(J2)=:WM(J2)+RP(l5L) 

WP(J1)=WP(J1)+RM(I,L) 

WM(J1)=WM(J1)+RP(I,L) 

862  CONTINUE 
3256  CONTINUE 

DO  897  J=1,NS 

897  WDOT(J,L)=(WP(J)“WM(J))/RHOOUT(L)/U(L) 

899  CONTINUE 
IOUT=IOUT+1 

63  IF(IFINIS)  65,69,65 
65  IF(X-XMAX)  67,66,66 

67  IF(PRNT-PCNT)  69,69,68 

68  CONTINUE 
GO  TO  5 
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66  IFINIS=2 
69  CALL  OUTPUT 
PCNT=0.0 

IF(IFINIS-I)  5,5,6 
C 

C  CHECK  DIFFUSION  STEP  SIZE 
C 

5  XD=DELPSI*DELPSI*SIGMA(1)/XMU(1)/XLE(1)/12.0  *FDL 
DO  511  I=2,NPSI 
DUMMY=A(I+1 )+A(I-1 )+A(I)+A(I) 

DUMMY=  PSI ( I ) *DELPSI*DELPSI *SIGMA ( I ) /XLE ( I ) /DUMMY/ 1 . 5*FDL 
511  XDrAMINKXD,  DUMMY) 

DX=AMIN1(DX,XD) 

DO  101  I=2,NPSI 

EX1  rPSKi)  »DELPSI**2/DX 

EX11=.5*(A(I)+A(I+1)) 

EX12=.5*(A(I)+A(I-1)) 

C 

C  INTEGRATE  MOMENTUM  EQUATION 
C 

RU(I)=(EX11*(U(I+1)-U(I))+EX12*(U(I-1)-U(I)))/EX1+U(I) 

RU(I)=RU(I)-DX*DPDX/RHO(I)/U(I) 

EX3=0.0 

EX4=0.0 

DO  21  J=1,NS 

EX3=EX3+H( J , I ) *WDOT( J ,  I ) 

21  EX4=EX4+CP( J,I)*(ALPHA( J, 1+1 )-ALPHA( J, 1-1 ) ) 

EX2rEX1*CPBAR(I) 

EX5=XLE(I)«A(I)/SIGMA(I) 

EX6=.5»(EX5+XLE(I+1 )*A(I+1 )/SIGMA(I+1 ) ) 

EX7=.5*(EX5+XLE(I-1 )«A(I-1 )/SIGMA(I-1 ) ) 

EXSrCPBAR ( I )*A( I )/SIGMA( I ) 

EX9= . 5*(EX8+CPBAR ( 1+1 ) *A( 1+1 )/SIGMA( 1+1 ) ) 

EX10=.5*(EX8+CPBAR(I-1 )*A(I-1 )/SIGMA(I-l ) ) 

EX14rEX4*EX5/4.0 

C 

C  INTEGRATE  ENERGY  EQUATION 

C 

RT( I ) = ( U( 1+1 ) -U( 1-1 ) ) ♦*2*A ( I )/EX2/ 4 .0+DX*DPDX/RHO( I )/CPBAR ( I ) +T( I ) 
1 +( ( EX9+EX1 4 ) *T ( 1+1 ) +( EX1 0-EX1 4 ) »T( 1-1 ) -(EX9+EX1 0 ) *T( I ) ) /EX2-EX3»DX 
2/CPBAR(I) 

RHOUIX=DX/(RHOOUT(I)*  U(I)) 

C 

C  INTEGRATE  SPECIES  EQUATIONS 

C 

DO  41  J=1,NS 

41  QXKJ)  =(EX6*(ALPHA(J,I+1)-ALPHA(J,I))+EX7»(ALPHA(J,I-1)-ALPHA 
1(J,I)))/EX1+ALPHA(J,I)+  QX(J,I)*RHOUIX 
DO  781  M=1 ,NS 
DO  781  N=1,NS 

CM1(M,N)=  CM(M,N,I)*RHOUIX 
IF  (M.EQ.N)  CM1(M,N)=CM1(M,N)  +1.0 
781  CONTINUE 

CALL  SLDP(QX1,CM1,NS) 

785  FORMAT  (1H  ,  215) 
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DO  782  J=1 ,  NS 

782  RALPHA(J,I)=  QXl(J) 

101  CONTINUE 

EX3=4.0«XMU( 1 )*DX/DELPSI/DELPSI 
RHOUIX=DX/(RHOOUT(1)*  U(1)) 

COMPUTE  U  AT  CENTER  LINE 

RU( 1 )=EX3*(U(2)-U( 1 ) )+U( 1 )-DX*DPDX/RHO( 1 )/U( 1 ) 

EX4=0.0 
DO  200  J=1 ,NS 
EX4rEX4+H( J , 1 ) *WDOT( J , 1 ) 

RALPHA( J, MPSI ) rALPHAC J, MPSI) 

200  QXKJ)  =EX3*XLE(1)*(ALPHA(J,2)-ALPHA(J,1))/SIGMA(1)+ALPHA(J,1) 
1+  QX(J,1)»RH0UIX 
DO  783  M=1,NS 
DO  783  N=1,NS 

CM1(M,N)=  CM(M,N,1)*RHOUIX 
IF  (M.EQ.N)  CM1(M,N)=CM1(M,N)  +1.0 

783  CONTINUE 
CALL  SLDP(QX1,CM1,NS) 

DO  784  J=1 ,  NS 

C  COMPUTE  SPECIES  AT  CENTER  LINE 
C 

784  RALPHA(J,1)=  QXI(J) 

C 

C  CALCULATE  TEMP.  AT  CENTER  LINE 
C 

RT( 1 ) =EX3*(T(2)-T( 1 ) )/SIGMA( 1 )+T( 1 )+DX»DPDX/RHO( 1 )/CPBAR( 1 ) 
1-EX4*DX/CPBAR(1) 

RT(MPSI)=T(MPSI) 

IF(IEDGE)  230,231,230 
C 

C  COMPUTE  TEMP.  AND  U  AT  EDGE 
C 

230  RU ( MPSI ) =U( MPSI ) -DX»DPDX/ RHO( MPSI ) /U( MPSI ) 

RT(MPSI) =T(MPSI)+DX»DPDX/RHO( MPSI ) /CPBAR(MPSI ) 

DO  210  I=MPSI,29 

RU(I)=RU(MPSI) 

U(I)=RU(MPSI) 

RT(I)=RT(MPSI) 

210  T(I)=RT(MPSI) 

231  CONTINUE 
1  IFINIS=1 

921  SAVEX=X 
SAVEDX=DX 
DO  941  1=1 ,29 
SAVEU(I)=U(I) 

SAVET(I)=T(I) 

DO  940  J=1 ,NS 
SAVEA(J,I)=ALPHA(J,I) 

940  CONTINUE 

941  CONTINUE 
MINIT  =13 
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MHALF  =  25 
NTEST=MPSI-1 
DO  967  I=1,NTEST 

CHECK  NEGATIVE  MDLE  FRACTION 

965  DO  967  J=1,NS 

IF  (RALPHA(J,I))  995,967,967 
967  CONTINUE 
X=X+DX 
PCNTrPCNT+DX 
DXrXD 

DO  925  1=1 ,29 
DO  926  J=1,NS 
926  ALPHA(J,I)=RALPHA(J,I) 

T(I)=RT(I) 

925  U(I)=RU(I) 

GO  TO  999 

995  IF  (DX.LT.DXMIN)  GO  TO  8000 

981  DX=SAVEDX/2.0 
X=SAVEX 

DO  985  1=1 ,29 
DO  982  J=1,NS 

982  ALPHA(J,I)=SAVEA(J,I) 

T(I)=SAVET(I) 

985  U(I)=SAVEU(I) 

GO  TO  2 

IF  MPSI  .GE.26  ,MPSI  IS  HALVED 
999  IF(MPSI-MHALF)  1001,1500,1500 

1001  IF(ABS(U(NPSI)-U(MPSI))/U(MPSI)-.010E0)  1011,1011,1004 
1011  IF(ABS(T(NPSI)-T(MPSI))/T(MPSI)-.050E0)  1002,1002,1004 

1002  CONTINUE 
GO  TO  2000 

1004  MPSI=MPSI+1 
NPSI=MPSI-1 
DO  1101  I=MPSI,29 
SAVEU(I)=U(NPSI) 

RU(I)=U(NPSI) 

U(I)=U(NPSI) 

SAVET(I)=T(NPSI) 

T(I)=T(NPSI)  ' 

RT(I)=T(NPSI) 

DO  1102  J=1,NS 
SAVEA(J,I)=ALPHA(J,NPSI) 

ALPHA( J, I)rALPHA( J,NPSI) 

1102  RALPHA(J,I)rALPHA(J,NPSI) 

1101  CONTINUE 
GO  TO  2000 
1500  IFINIS=0 

DELPSI=DELPSI+DELPSI 
DO  1600  I=1,MINIT 
DO  1650  J=1,NS 

1650  ALPHA(J,I)=ALPHA(J, 2*1-1) 
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T(I)=T(2»I-1) 

1600  U(I)=U(2«I-1) 

MPSI=MINIT 
NPSI=MPSI-1 
DO  1700  I=MINIT,29 
DO  1750  J=1,NS 
ALPHA(J,I)=ALPHA(J,MPSI) 

1750  RALPHA(J,I)=ALPHA(J,MPSI) 

T{I)=T(MPSI) 

RT(I)=T(MPSI) 

U(I)=U(MPSI) 

1700  RU(I)=U(MPSI) 

DO  1800  1=2,29 
1800  PSI(I)=PSI(I-1)+DELPSI 
ITER=0 
ISTEP=0 
GO  TO  2000 

8000  WRITE(NDBG,8001) 

8001  FORMAT(68H1 NEGATIVE  PARAMETER  -  NOT  CORRECTED  BY  REPEATED  HALVING 
10F  STEP  SIZE) 

IFINIS=2 
GO  TO  69 
2000  CONTINUE 

CALL  TICK(ISECS) 
lELAPS  =  ISECS-ISECST 
IF(IELAPS,LT.0)  lELAPS  =  IDIFFT  +  ISECS 
IF(IELAPS.GE.ILIMIT)  GO  TO  6 
GO  TO  2 

100  FORMAT(14I5) 

102  FORMAT(8F10.4) 

111  FORMAT  (7(  E10.3)) 

222  FORMAT(A6,7E10.3) 

333  FORMAT(10A8) 

444  F0RMAT(A6,1X,A6,8X,A6,1X,A6,1X,A6,7X,I2,I1,E8.2,F4.1,F9.1) 

555  F0RMAT(8(  El 0.3)) 

666  FORMAT(10I5) 

1000  FORMAT(7E10.3) 

9900  FORMAT  (39H1  MIXING  REGION  INTERSECTS  AXIS  AT  X  =  El 5.7) 

6  CONTINUE 
XORJ=X/RJ 

790  FORMATC  E10.3,60X,E10.3) 

IF  (IPUNCH  .EQ.  0  .OR.  IPUNCH  .EQ.  1)  GO  TO  9 
X=X«AMULT 
XMAX=XMAX^AMULT 
PRNT=PRNT*AMULT 
RJ=RJ«AMULT 
DXMIN=DXMIN«AMULT 
PC(2)=PC(2)/AMULT 
PC(3)=PC(3)/AMULT1 
PC(4)=PC(4)/AMULT2 
DELPSI=DELPSI*AMULT3 
DO  7666  I=1,MPSI 
7666  U(I)=U(I)*AMULT 

WRITE ( NUNITE, 333 ) ( TITLE( I ) , 1= 1 , 1 0) 

WRITE ( NUNITE, 666 )MPSI, NS, ITURB, NR, IOUT1 , I0UT2, IPUNCH, ITIME, 
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*IPRESS,NT 

WRITE(NUNITE,1000)FRE]QA(1 )  ,FREQA(2)  ,FRBQA(3)  ,FRBQA(i|)  ,FRE1QA(5) , 
*FREQA(6) 

WRITEC NUNITE, 1000)  X, XMAX, PRNT, XLE( 1 ) , SIGMA( 1 ) , RJ , XK2 
WRITE(NUNITE,1000)  DXMIN,  FDL,  PC(1 ),PC(2),PC(3),PC(4) 

WRITEC  NUNITE,  1000)  PPUNCH,T(  1 )  ,T(MPSI)  ,U(  1 )  ,U(MPSI)  ,DELPSI,  THNET 
WRITEC NUNITE, 1000)  CTCD , 1=1 ,MPSI) 

WRITE C NUNITE, 1000)  C UC I ) , 1= 1 , MPSI ) 

DO  8  I  =  1,MPSI 
DO  1104  J=1,NS 

1104  RALPHACJ,I)=RALPHACJ,I)/WTMIXCI) 

WRITE C  NUNITE , 1 000 )  C RALPH AC  J , I ) , J= 1 , NS ) 

8  CONTINUE 

9  CONTINUE 
STOP 
END 

SUBROUTINE  OUTPUT 

DIMENSION  AC30) ,RHOC30) ,YC30) ,TC30) ,PSIC30) ,RTC30) ,SUMC30) ,ARC25) , 
1HSTATC30) , HC25 , 30) , ALPHAC25 , 30) , RALPHAC25 ,30) , CPC25 , 30) , SIGMAC30) , 
2  WIMOLEC25) ,CPBARC30) ,CC25,9) ,AIDC25) ,ETAC30) ,RATIOC30) , 

3RUC30),UC30),TITLEC12),XLEC30),XMUC30)  ,GC25)  ,WTMIXC30) , 

4RCC49,3),IRRRC49,5),WPC25),WMC25),WDOTC25,30),SAVETC30),SAVEU(30), 

5  IRRC49),FREQC30),SAVEAC25,30),  PCC4),ZIDC5) , 

6  ECCC30) ,HOUTC30) ,yOUTC30) ,RHOOUTC30) ,XMUOUTC30) ,XLTC30) , 
7T4C30),TFDGC30),IRTC49),RPC49,30),RMC49,30) 

DIMENSION  ISAVEC6),  FREQAC6),  ALOCC50,6),  ATTC6),  YATTC50) 
DIMENSION  ZCONC 16,30) 

COMMON/CONSTS/AMULT,AMULT1 ,AMULT2,AMULT3 
COMMON/UNITS/NUNITA, NUNITB, NUNITC, NUNITD, NUNITE, NUNITF, 

*NUNITG, NUNITH , NUNITI, NUNITJ , NUNITK, NUNITL, NUNITM, NUNITN , 
*NOUT,NDBG,NNNOUT 
COMMON/C/  IZSPEC,ISPECC16) 


COMMON 

A  , 

RHO 

f 

Y 

9 

T  , 

PSI  , 

RT 

COMMON 

SUM  , 

AR 

HSTAT 

9 

H  , 

ALPHA  , 

RALPHA 

COMMON 

CP  , 

SIGMA 

9 

WTMOLE  , 

CPBAR  , 

C 

COMMON 

AID  , 

ETA 

1 

RATIO 

9 

RU  , 

u  , 

TITLE 

COMMON 

XLE  , 

XMU 

9 

G  , 

WTMIX  , 

WDOT 

COMMON 

SAVEU  , 

SAVET 

1 

WM 

9 

WP  , 

RC 

COMMON 

SAVEA 

9 

PC  , 

X  , 

XMAX 

COMMON 

PRNT  , 

DXMIN 

9 

DX 

9 

FDL  , 

DELPSI  , 

RJ 

COMMON 

XK2  , 

P 

9 

ZID 

9 

FREQ  , 

ECC  , 

DPDX 

COMMON 

Y  OUT 

,  HOUT 

,  RHOOUT  ,  IRRR 

,  IRR 

,  ifin; 

COMMON 

IPAGE  , 

MPSI 

9 

MY 

9 

NS  , 

NR  , 

lEDGE 

COMMON 

ITURB  , 

IPRESS 

9 

NPSI 

9 

ITEST  , 

ITER  , 

lECC 

COMMON 

IRT  , 

XMU 

OUT  ,  XLT 

,  T4 

,  TFDG 

,  lOUT 

COMMON  IOUT1  ,  I0UT2  ,RP  ,  RM  ,  ISAVE  ,  IPUNCH 
COMMON  TKINET, NFREQA, ALOC, FREQA, QQ1 00 , QQ200 , QQ300 , QQ400 
DATA  AMULT5/1 000 . / , AMULT4/1 . 488/ 

DATA  BLANK/8H  / 

DATA  ZCO  /6HCO  / 

DATA  ZC02/6HC02  / 

DATA  ZH20/6HH20  / 

DATA  ZO  /6HO  / 

DATA  ZOH/6HOH  / 

DATA  ZAL/6HAL203  / 
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DUMMY  =  0.0 

BIG=1.0E30 

NS1=NS+1 

IF(NS1  .GT.  25)  GO  TO  2 
DO  1  I=NS1,25 

1  AID(I)=BLMK 

2  DO  5  1=1 ,  NS 

IF  (AID(I).BQ.ZCO  )  ICO  =I 
IF  (AID(I).EQ.ZC02)  IC02=I 
IF  (AID(I).EQ.ZH20)  IH20=I 
IF  (AID(I).EQ.ZO  )  10  =I 

IF  (AID(I).BQ.ZOH)  IOH=I 
IF  (AID(I).EQ.ZAL)  IAL=I 
5  CONTINUE 
IF(IECC)  531,539,531 

531  DO  532  I=1,MPSI 

532  ECC( I ) =RHO ( I ) *ALPHA( lECC, I) *3 . 1 08E23 
539  DO  10  I=1,MPSI 

yOUT(I)=Y(I)/RJ 
XMU0UT(I)=XMU(I)*32.174 
HOUT( I ) =HSTAT( I )/45055 .31 
SUM(I)=0.0 
DO  10  J=1,NS 

10  SUM(I)=SUM(I)+ALPHA(J,I)*WTMOLE(J) 

UD=.05*U(1)+.95*U(MPSI) 

DO  83  I=2,MPSI 

IF  ((U(I)-UD)*(U(I-1)-UD))  84,84,83 

83  CONTINUE 

84  VR=(Y(I)-Y(I-1 ) )*(UD-U(I-1 ) )/(U(I)-U(I-1 ))+Y(I-1 ) 
TD=.05»T(1)+.95*T(MPSI) 

DO  85  I=2,MPSI 

IF  ((T(I)-TD)*(T(I-1)-TD))  86,86,85 

85  CONTINUE 

86  TR=(Y(I)-Y(I-1 ) )*(TD-T(I-1 ) )/(T(I)-T(I-1 ) )+Y(I-1 ) 

TR=TR/RJ 

VRrVR/RJ 
DO  87  J=1,NS 
AR(J)=0.0 

AD= .05*ALPHA( J , 1 ) + .95*ALPHA( J , MPSI ) 

IF(ALPHA(J,MPSI))  91,92,91 

91  DO  88  I=1,MPSI 
IF(ALPHA(J,I)-AD)  88,88,89 

88  CONTINUE 

89  AR( J)=Y(I-1 )+(Y(I)-Y(I-1 ) )*(AD-ALPHA( J, 1-1 ) )/(ALPHA( J, I)-ALPHA( J, I 
1-1)) 

AR(J)=AR(J)/RJ 
GO  TO  87 

92  DO  93  I=1,MPSI 
IF(ALPHA(J,I)-AD)  94,93,93 

93  CONTINUE 

94  GO  TO  89 

87  CONTINUE 
PCNT=0.0 
IPAGE=IPAGE+1 
XXOUT=X*AMULT 
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XXX=XXOUT 

WRITE  (NOUT,201)XXOUT,(TITLE(I),I=1,10),IPAGE 
WRITE  (NOUT,102) 

XORJ=X/RJ 

POUT=P/2117.0 

DPOUT=DPDX/2117.0 

DDX=DX«AMULT 

WRITE(NOUT,103)  XORJ,DDX,POUT 
IF  (ITURB-6)  8600,8500,8600 
8500  WRITE  (N0UT,8555) 

QQ101=QQ100/RJ 

QQ201=QQ200/RJ 

WRITE  (NCXJT.8556)  QQ101  ,QQ201  ,QQ300,C3Q400 

8555  FORMAT(1H0,8X,4HHALF,21X,12HINNER  MIXING, 17X, 11 HMACH  NUMBER, 16X,11 
*HMIXING  RATE/7X,10HRADIUS/R  ,16X,15HZ0NE  RADIUS/ R  ,14X,14HAT  HAL 
*F  RADIUS, 1 4X, 1 1 HCOEFFICIENT) 

8556  FORMAT  (4X,  E14.6,  3E28.6) 

8600  WRITE  (NOUT,107) 

WRITE  (NOUT,509) 

TTT=T(1) 

WRITE(NNNOUT, 8888)XXX, TTT 
8888  FORMAT(F6.1,F8.1) 

DO  73  I=1,MPSI 

SS1=  89517. 501*WTMIX(I) 

SS2=  CPBAR(I)/(CPBAR(I)-SS1) 

SSrSQRT(SS2*SS1*T(I)) 

XMACH=  U(I)/SS 
UU=U(I)«AMULT 
RRHOUT=RHOOUT(I ) *AMULT5 
XXM0UT=XMU0UT( I ) *AMULT4 
PPSI=PSI(I)*AMULT3 
IF(IECC)  71,72,71 

71  WRITE  (NOUT,207)I,YOUT(I),UU,T(I),RRHOUT,XMACH,  HOUT(I), 
«XXMOUT,ECC(I),PPSI,I 

GO  TO  73 

72  WRITE  (NOUT,307)I,YOUT(I),UU,T(I),RRHOUT,XMACH,  HOUT(I), 
*XXMOUT,PPSI,I 

73  CONTINUE 

DO  581  I=1,MPSI 
DO  581  J=1,NS 

581  RALPHA(J,I)=ALPHA(J,I)/WTMIX(I) 

IRPT=(NS+6)/7 
DO  564  KK=1,IRPT 
I1=1+(KK-1)*7 
I2=7+(KK-1)*7 

WRITE  (NOUT,201 )XXOUT,(TITLE(I),I=1 ,10),IPAGE 
WRITE  (NOUT,409) 

IF  (I2.GE.25)  GO  TO  50 

WRITE  (NOUT,108)(AID(J),J=I1,I2) 

DO  81  I=1,MPSI 

81  WRITE  (NOUT,208)I,YOUT(I),(RALPHA(J,I),J=I1,I2),I 
IF(IOUT1)564,564,74 

74  WRITE  (NOUT,420) 

WRITE  (N0UT,421)(AID(J),J=I1,I2) 

DO  82  I=1,MPSI 
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IF(T(I)-TKINET)  564,564,82 
82  WRITE  (N0UT,422)I,(WD0T(J,I),J=I1,I2),I 

GO  TO  564 
50  12=25 

WRITE  (NOUT,108)  (AID(J) , J=I1 ,12) 

DO  52  I=1,MPSI 

52  WRITE  (NOUT,53)  I,Y0UT(I),(RALPHA(J,I),J=I1,I2),I 

53  FORMAT  (I3,F9.5,4E13.5,42X,I3) 

IF  (IOUT1)  564,564,54 

54  WRITE  (NOUT,420) 

WRITE  (NOUT,55)  (AID(J) , J=I1 ,12) 

55  FORMAT  (3H0PT,8X,4(3X,A6,4X),43X,3H  PT) 

DO  56  I=1,MPSI 

IF  (T(I)-TKINET)  564,564,56 

56  WRITE  (NOUT,57)  I,(WDOT(J,I),J=I1 ,I2),I 

57  FORMAT  (I3,9X,4E13.5,42X, 13) 

564  CONTINUE 

IF(IOUT2)567,567,75 
75  IRPT=(NR+9)/10 
N=0 

NNR=NR-1 

DO  565  KK=1,IRPT 

LL=0 

N=N+1 

WRITE  (NOUT,201)XXOUT,(TITLE(I),I=1,10),IPAGE 

65  I1=1+(N-1)«5 
I2=5+(N-1)*5 
NNN1=I1 
NNN2=I1+1 
NNN3=I1+2 
NNN4=I1+3 
NNN5=I2 

WRITE(NOUT,209) 

WRITE  ( NOUT, 43 1 ) NNN1 , NNN2 , NNN3 , NNN4 , NNN5 
WRITE  (NOUT, 432) 

DO  63  I=1,MPSI 
IF(T(I)-TKINET)  566,566,63 

63  WRITECNOUT, 433)1, YOUT(I) , (RP( J, I) ,RM( J,I) , J=I1 ,12) ,I 

566  IF(NNR/(5*N))565,565,64 

64  IF(LL)565,66,565 

66  N=N+1 
LL=1 

GO  TO  65 
565  CONTINUE 

567  CONTINUE 
568  CONTINUE 

1065  CONTINUE 

WRITE  (NOUT, 1068) 

1066  DO  602  I=1,NPSI 

FT1  =  1.0/SQRT(T(I)) 

FT2  =  1.0/FT1 
FT3  =  1.0/T(I) 

FT4  =  T(I)*»0.75 
SUMS  =0.0 
DO  603  IDX  =  1,6 
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IF(ISAVE(IDX).EQ.0)  GO  TO  603 
K  =  ISAVE(IDX) 

TERM  =  RALPHA(K,I) 

GO  TO  (604, 605, 606, 607, 608, 614), IDX 

604  Q  =  (1.29E-17)*FT2  +  2.46E-16 
GO  TO  609 

605  Q  =  (0.758E-13)*Fn 
GO  TO  609 

606  Q  =  (1.53E-11)*FT3 
GO  TO  609 

607  Q  =  (9.0E-18)*FT2  +  8.9E-16 
GO  TO  609 

608  Q  =  3.29E^23*6.21E5*FT2 
GO  TO  609 

614  Q=  1.85»(6.21)**(-2)»(1.0E-10)*FT3 

609  SUMS  =  SUMS  +  Q*TERM 
603  CONTINUE 

XNEU  =  (4.57E27)*SUMS*POUT»FT1 
ECON=  0.07157*  ECC(I)/XNEU/.0254 

602  IF(IZSPEC.EQ.0)  GO  TO  613 
WRITE(NUNITD,786)  X0RJ,MPSI 
DO  615  I=1,MPSI 
Yy=y(i)/Rj 
PP=P/2117.0 

WRITE(NUNITD,785)  YY,T(I),PP 
DO  616  M=1,nSPEC 
K=ISPEC(M) 

616  ZCON(M,I)=RALPHA(K,I) 

LINES±IZSPEC/7+1 .1 
IF(IZSPEC.EQ.7)  LINESrI 
DO  559  L=1, LINES 
LSTART=(U1  )*7+1 
LENDrMIN0  ( IZSPEC,  L*7) 

559  WRITE(NUNITD,785)  (ZC0N(K,I) ,K=LSTART,LEND) 

615  CONTINUE 
613  CONTINUE 

102  FORMAT(/,8X,3HX/R,8X,14HDELTA  X  METERS, 4X,10HPRESS( ATM)) 

103  FORMAT(4X,  6E14.6) 

107  FORMAT(4H0  PT,5X,3HY/R.  6X,8HVEL0CITY,4X,11HTEMPERATURE,5X,7HDENSI 
1TY,  4X8HMACH  NO.  ,  8X,8HEN'fflALPY,5X,9HVISCOSITY, 

2  8X,3HPSI,12X,2HPT  ) 

108  FORMAT(103X,3H  PT,T1 ,3H0PT,3X,5H  Y/R  ,7(3X,A8,2X)) 

201  FORMAT(1H1,//////3H  X=  E15.7,7H  METERS, 8X,10A8,8X,4HPAGEI4) 

207  FORMAT(I4,F10.4,  8E14.6,I4) 

208  FORMAT(I3,F9.5,  7E13.5,I3) 

209  FORMATCIH  ,//40X,28HREACTION  RATES  (MOLE/ML-SEC)//) 

307  FORMAT(I4,F10.4,  7E14.6,4X,I4) 

409  FORMAT(1H0,44X,14HMOLE  FRACTIONS) 

420  FORMAT(1H0,35X,36HNET  RATE  OF  PRODUCTION  (W-DOT/RHO*U)) 

421  FORMAT(3H0PT,8X,7(3X,A6,4X),1X,3H  PT) 

422  F0RMAT(I3,9X,  7E13.5,I3) 

431  FORMAT(1H0,2HPT,4X,3HY/R,  8X,5(8HREACTI0N,I3,11X),2HPT) 

432  F0RMAT(19X,5(2HRP,  9X,2HRM,10X)) 

433  FORMAT(I3,1X,  11E11.4,I4) 

509  FORMAT(18X,10HMETERS/SEC,4X,10H  K  ,6X,5HGM/CC,22X,6HCAL/GM, 
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*7X,8HKG/M/SEC) 

610  F0RMAT(3X,I3,5X,F7.3,8(2X,  E10.4)) 

611  F0RMAT(3E15.7) 

612  F0RMAT(E15.7,I5) 

780  F0RMAT(8(  El 0.3)) 

785  FORMAT(10X,7E10.3) 

786  FORMAT(E10.3,I10) 

790  FORMATC  E10.3,60X,E10.3) 

1068  FORMAT  (1H0) 

RETURN 

END 


SUBROUTINE  INOUT  ^  ^ 

DIMENSION  A(30),RHO(30),Y(30),T(30),PSI(30),RT(30),SUM(30),AR(25); 
1HSTAT(30),H(25,30),ALPHA(25,30),RALPHA(25,30),CP(25,30),SIGMA(30); 
2  WTMOLE(25),CPBAR(30),C(25,9),AID(25),ETA(30),RATIO(30), 

3RU(30) ,U(30) ,TITLE(12) ,XLE(30) , XMU(30)  ,G(25) ,WTMIX(30) : 

4RC(49,3) ,IRRR(49,5) ,WP(25) ,WM(25) ,WDOT(25,30) ,SAVET(30) ,SAVEU(30)  i 

5  IRR(49),FREQ(30),SAVEA(25,30)i  PC(4),ZID(5), 

6  ECC(30) ,HOUT(30) ,YOUT(30) ,RHOOUT(30) ,XMUOUT(30) ,XLT(30)  ■ 
7T4(30),TFDG(30),IRT(49),RP(49,30),RM(49,30) 

DIMENSION  ISAVE(6) 


COMMON/CONSTS/ AMJLT , AMULT1 , AMULT2 , AMULT3 
COMMON/UNITS/  NNDUM( 14), ND, NDBG 


COMMON 

A 

RHO 

^  # 

,  Y 

,  T 

> 

PSI 

,  RT 

COMMON 

SUM 

t 

AR 

,  HSTAT 

r  H 

ALPHA 

,  RALPHA 

COMMON 

CP 

1 

SIGMA 

,  WTMOLE 

j 

CPBAR 

,  c 

COMMON 

AID 

1 

ETA 

,  RATIO 

,  RU 

1 

U 

,  TITLE 

COMMON 

XLE 

f 

XMU 

> 

G 

9 

WTMIX 

,  WDOT 

COMMON 

SAVEU 

f 

SAVET 

,  WM 

,  WP 

9 

RC 

COMMON 

SAVEA 

t 

PC 

9 

X 

,  XMAX 

COMMON 

PRNT 

} 

DXMIN 

,  DX 

,  FDL 

9 

DELPSI 

,  RJ 

COMMON 

XK2 

P 

,  ZID 

,  FREQ 

9 

ECC 

,  DPDX 

COMMON 

YOUT 

f 

HOUT 

,  RHOOUT 

,  IRRR 

9 

IRR 

,  IF  INIS 

COMMON 

IPAGE 

f 

MPSI 

,  MY 

,  NS 

9 

NR 

,  lEDGE 

COMMON 

ITURB 

f 

IPRESS 

,  NPSI 

,  ITEST 

9 

ITER 

,  lECC 

COMMON 

IRT 

f 

XMUOUT 

,  XLT 

,  T4 

9 

TFDG 

,  lOUT 

COMMON 

I0UT1 

f 

I0UT2 

,RP 

,  RM 

9 

ISAVE 

,  IPUNCH 

COMMON  TKINET 


X=X«AMULT 
XMAXrXMAX*AMULT 
PRNT=PRNT«AMULT 
RJ=RJ*AMULT 
DXMINrDXMIN»AMULT 
DO  7666  I=1,MPSI 
7666  U(I)=U(I)*AMULT 
C 

WRITE(ND,1901) 

WRITE(ND,1902) 

WRITE ( ND, 1 903 ) ( TITLE( I ) , 1= 1 , 1 0 ) 

IF(IPRESS.NE.0)  WRITE(ND,1904)P 

IF(IPRESS.EQ.0)  WRITE(ND,1905)P 

WRITE(ND,1908)RJ 

WRITE( ND, 1 909)XLE( 1 ) , SIGMA( 1 ) 

WRITE(ND,1966)X,XMAX 
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WRITE(ND, 1 967) PRNT, DXMIN 
LL=ITURB+2 

GO  TO  (1991, 1999, 1945, 1978, 1926, 1933, 9950, 8000), LL 
1991  WRITE(ND,1959) 

GO  TO  2999 

1999  WRITE(ND,1960)XK2 
GO  TO  2999 

1945  WRITE(ND,1961)XK2 
GO  TO  2999 

1978  WRITE(ND,1962)XK2 
GO  TO  2999 
1926  GO  TO  2999 
1933  GO  TO  2999 

8000  WRITE(ND,8001) 

8001  FORMAT(1H0,22X,30HDONALDSON/GRAY  VISCOSITY  MODEL) 

GO  TO  2999 

9950  WRITE(ND,9951)XK2 
2999  CONTINUE 

WRITE(ND,1916) 

WRITE(ND, 1 906)T( 1 ) , T(MPSI) 

WRITE ( ND, 1 907) U ( 1 ), U ( MPSI ) 

WTMIX(1)=0.0 

WTMIX(MPSI)=0.0 

DO  1930  J=1,NS 

WTMIXd  )=WTMIX(1  )+ALPHA(J,1 ) 

1 930  WTMIX(MPSI)=WTMIX(MPSI)+ALPHA(J, MPSI) 

DO  1919  J=1,NS 

RALPHA(J,1  )=ALPHA(J,1)AmiIX(1) 

RALPHAC J , MPSI)=ALPHA( J , MPSI)/WTMIX(MPSI) 

1919  WRITE(ND,1917)AID(J),RALPHA(J,1),RALPHA(J,MPSI) 

WRITE(ND,120) 

DO  159  1=1, NR 
LrIRRd) 

GO  TO(131, 132, 133, 134, 135, 136, 137, 138, 139, 140), L 

131  J1=IRRR(I,1) 

J2=IRRR(I,2) 

J3=IRRR(I,3) 

J4=IRRR(I,4) 

WRITE(ND,121)I,AID(J1),AID(J2),AID(J3),AID(J4),(RC(I,J),J=1,3) 

GO  TO  159 

132  J1=IRRR(I,1) 

J2=IRRR(I,2) 

J3=IRRR(I,3) 

WRITE(ND,122)I,AID(J1),AID(J2),AID(J3),(RC(I,J),J=1,3) 

GO  TO  159 

133  J1=IRRR(I,1) 

J2=IRRR(I,2) 

J3=IRRR(I,3) 

J4rIRRR(I,4) 

J5=IRRR(I,5) 

WRITECND, 123)1, AID(J1 ) ,AID(J2) ,AID(J3) ,AID(J4) ,AID(J5) , (RC(I, J) , J 
11,3) 

GO  TO  159 

134  J1=IRRR(I,1) 

J2=IRRR(I,2) 
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J3=IRRR(I,3) 

WRITE(ND,124)I,AID(J1),AID(J2),AID(J3),(RC(I,J),J=1,3) 

GO  TO  159 

135  J1=IRRR(I,1) 

J2=IRRR(I,3) 

J3=IRRR(I,4) 

WRITE(ND, 125)1, AID(J1 ) ,AID(J2) ,AID(J3) ,(RC(I,  J) , J=1 ,3) 

GO  TO  159 

136  J1=IRRR(I,1) 

J2=IRRR(I,2) 

J3=IRRR(I,3) 

J4=IRRR(I,4) 

WRITECND, 126)1, AID(J1 ) ,AID( J2) , AID( J3) , AID( J4) , (RC(I, J) , J=1 ,3) 

GO  TO  159 

137  J1=IRRR(I,1) 

J2=IRRR(I,2) 

J3=IRRR(I,3) 

WRITECND, 127)1, AID(J1 ) ,AID( J2) , AID( J3) , (RC(I, J) , J=1 ,3) 

GO  TO  159 

138  J1=IRRR(I,1) 

J2=IRRR(I,2) 

J3=IRRR(I,3) 

J4=IRRR(I,4) 

J5=IRRR(I,5) 

WRITECND, 128)1, AID(J1 ) ,AID( J2) ,AID(J3) ,AID(J4) ,AID(J5) ,(RC(I, J) , J= 
11,3) 

GO  TO  159 

139  J1=IRRR(I,1) 

J2=IRRR(I,2) 

J3=IRRR(I,3) 

WRITECND, 129)1, AID(J1 ) ,AIDC J2) ,AIDCJ3) , CRCCI, J) , J=1 ,3) 

GO  TO  159 

140  J1=IRRRCI,1) 

J2=IRRRCI,2) 

J3=IRRRCI,3) 

WRITECND, 130)1, AIDCJ1),AIDCJ2),AIDCJ3), CRCCI, J),J=1, 3) 

159  CONTINUE 

120  FORMATC1H0,19X,26HREACTIONS  BEING  CONSIDERED, 6X, 19HKR=A*EXPCB/RT)/ 
1T»*N,7X,1HA,8X,1HN,9X,1HB,7X,23HCMOLECULE^ML-SEC  UNITS)) 

121  FORMATCI9,9X,A6,2H+  ,A6,8X,2H=  ,A6,2H+  ,A6,l8X,1E10.4,2X 
1,F4.1,2X,F10.1) 

122  FORMATCI9,9X,A6,2H+  ,A6,3H+  M,5X,2H=  ,A6,3H+  M,23X, 
1E10.4,2X,F4.1,2X,F10,1) 

123  FORMATCI9,9X,A6,2H+  ,A6,8X,2H=  ,A6,2H+  ,A6,2H+  ,A6,10X, 
1E9.3,2X,F4J,2X,F10.1) 

124  FORMATCI9,9X,A6,2H+  ,A6,8X,2H=  ,A6,26X,E9.3,2X,F4.1 ,2X,F10.1 ) 

125  FORMAT C 19, 9X,A6,3H+  M,13X,2H=  ,A6,2H+  ,A6,3H+  M,15X, 
1E9.3,2X,F4.1,2X,F10.1) 

126  FORMATCI9,9X,A6,2H+  ,A6,8X,2H=  ,A6,2H+  ,A6,18X,E9.3,2X, 

1 F4 . 1 , 2X, FI 0 . 1 , 3X , 1 6HONE  WAY  REACTION) 

127  FORMATCI9,9X,A6,2H+  ,A6,3H+  M,5X,2H=  ,A6,3H&  M,23X, 

1 E9 . 3 , 2X, F4 . 1 , 2X, FI 0 . 1 , 3X, 1 6HONE  WAY  REACTION) 

128  FORMATC  I9,9X,A6,2H+  ,A6,8X,2H=  ,A6,2H+  ,A6,2H+  ,A6,10X, 

1  E9.3,2X,F4.1,2X,F10.1,3X,16HONE  WAY  REACTION) 

129  F0RMATCI9,9X,A6,2H+  ,A6,8X,2Hr  ,A6,26X,E9.3,2X,F4.1 ,2X,F10.1 ,3 
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1X,16H0NE  WAY  REACTION) 

130  FORMAT(I9,9X,A6,3H+  M,13X,2H=  ,A6,2H+  ,A6,3H+  M,15X,E9.3, 
12X,F4.1,2X,F10.1,3X,16HONE  WAY  REACTION) 

1901  P0RMAT( 1 HI, 37X,46HAEROCHEM  RESEARCH  LABORATORIES  PRINCETON  N.J.) 

1902  FORMAT(35X,50HAXISYMMETRIC  MIXING  WITH  NON-EQUILIBRIUM  CHEMISTRY) 

1 903  FORMAT( 1 H0 , 24X, 1 0A8 ) 

1904  FORMAT(1H0,22X,19HPRESSURE(INITIAL)  =  E15.7,12H  ATMOSPHERES) 

1905  FORMAT(1H0,22X,20HPRESSURE(CONSTANT)  =  E15.7,12H  ATMOSPHERES) 

1906  FORMAT(23X,24HTEMPERATURE(DEG.  KELVIN), 3X,  E15.7,4X,  E15.7) 

1907  FORMAT(23X,24HVELOCITY  (METERS/SECOND) ,3X,  E15.7,4X,  E15.7) 

1908  FORMAT(1H0,22X,14HNOZZLE  RADIUS=  E15.7,7H  METERS) 

1909  FORMAT(1H0,22X,23HLEWIS  NUMBER ( CONSTANT )=  E15.7,5X,25HPRANDTL  NUM 
1BER( CONSTANT) =  El 5. 7) 

1 91 6  FORMAK 1 H0 , 54X, 3HJET, 1 6X, 4HEDGE) 

1917  FORMAT(23X,13HMOLE  FRACTION, 3X,A6,5X,  E15.7,4X,  E15.7) 

1959  FORMAT(1H0,22X,40HUMINAR  VISCOSITY  M0DEL( SUTHERLANDS  LAW)) 

1960  FORMAT(1H0,22X,29HCONSTANT  VISCOSITY  MODEL  MU=  El 5. 7) 

1962  FORMAT(1H0,22X,31HTING-LIBBY  VISCOSITY  MODEL  K=  El 5. 7) 

1961  FORMAT(1H0,22X,27HFERRI  VISCOSITY  MODEL  K=  E15.7) 

1966  F0RMAT(/,22X,18HX  INITIAL(METERS)=,E15.7,12X,16HX  FINAL( METERS) =, 
*E15.7) 

1967  FORMAT(1H0,22X,16HPRINT  INCREMENT=  E15.7,12X,18HMINIMUM  STEP  SIZE 
1=  E15.7) 

9951  FORMAT(1H0,22X,69HTING-LIBBY  VISCOSIH  MODEL  AFTER  MIHNG  REGION  I 
INTERSECTS  X  AXIS  K=  El 5. 7) 

X=X/AMULT 
XMAXrXMAX/AMULT 
PRNT=PRNT/AMULT 
RJ=RJ/AMULT 
DXMIN=DXMIN/AMULT 
DO  7667  I=1,MPSI 
7667  U(I)=U(I)/AMULT 
RETURN 
END 

SUBROUTINE  GRATE(ANSWER, Y,X,N) 

DIMENSION  X(30),Y(50) 

COMMON/UNITS/NNDUM( 1 5 ) , ND 
SUM  =  0.0 
INTER  =  N-1 

11  =  1 

12  =  2 

DO  1  I  z  1, INTER 

SUM  =  SUM  +  (X(I2)-X(I1))»(Y(I2)+Y(11)) 

11  =  II  +  1 

12  =  12  +  1 

1  CONTINUE 

ANSWER  z  0.5»SUM 

RETURN 

END 

SUBROUTINE  TICK(JJJJ) 

CALL  SECOND(TIME) 

JJJJzTIME 

RETURN 

END 

SUBROUTINE  SLDP(X,A,N) 
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THIS  PROGRAM  FINDS  THE  SOLUTIONS  TO  A  SET  OF  N  SIMULTANEOUS  LINEAR 
EQUATIONS  BY  USING  THE  GAUSS-GORDAN  REDUCTION  ALGORITHM  WITH  THE 
DIAGONAL  PIVOT  STRATEGY 
DIMENSION  A(25,25),X(25) 

COMMON/UNITS/  NNDUM( 1 4 ) , ND, NDBG 
DO  9  K=1,N 

IF  (ABS(A(K,K))  .GT.  1.E-10)  GO  TO  5 
WRITE(NDBG,101) 

GO  TO  99 

5  KP1=  K+1 

DO  6  J=  KP1,  N 

6  A(K,J)=  A(K,J)/A(K,K) 

X(K)=  X(K)/A(K,K) 

A(K,K)=  1.0 

DO  9  1=1, N 

IF  (I.EQ.K  .OR.  A(I,K).EQ.0.)  GO  TO  9 
DO  8  J=KP1,N 

8  A(I,J)=  A(I,J)-  A(I,K)*A(K,J) 

X(I)=  X(I)~  A(I,K)*X(K) 

A(I,K)=0. 

9  CONTINUE 
99  CONTINUE 

101  FORMATC  22H  ERROR -  SMALL  PIVOT  ) 

RETURN 

END 

SUBROUTINE  TKEY(T, TIB, ITKEY, SDT, HDT, NT, J) 

DIMENSION  TTB(30,24) 

COMMON/UNITS/NNDUM( 15), NDBG 

NT1=NT-.1 

DO  10  IT=1,NT1 

DT=  TTB(IT+1,J)-TTB(IT,J) 

SDT=(T-TTB(IT,J))/DT 
HDT=(TTB(  IT+1,J)-T)  /DT 
IF  ((SDT*HDT).GE.0.0)  GO  TO  11 

10  CONTINUE 
WRITE(NDBG,100)  T,IT 
ITKEY=0 

100  FORMATdH  ,  28H  TEMPERATURE  OUT  OF  RANGE  ,  El 4. 5, 15) 

RETURN 

11  ITKEY=IT 
RETURN 
END 

SUBROUTINE  LIPLN( ITKEY, I, ATB, SDT, HDT, AX) 

DIMENSION  ATB(25,30) 

AX=  ATB(I,ITKEY)»HDT+  ATB(I, ITKEY+1 )»SDT 

RETURN 

END 
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APPENDIX  K 

MEFF,  BLAKE,  LAPP  OUTPUT  LISTING 
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1 .  BRL  Report  Number _ ^Date  of  Report 


2.  Date  Report  Received 


3.  Does  this  report  satisfy  a  need?  (Comment  on  purpose,  related  project,  or 
other  area  of  interest  for  which  the  report  will  be  used.) 


4.  How  specifically,  is  the  report  being  used?  (Information  source,  design 
data,  procedure,  source  of  ideas,  etc. ) _ 


5.  Has  the  information  in  this  report  led  to  any  quantitative  savings  as  far 
as  man-hours  or  dollars  saved,  operating  costs  avoided  or  efficiencies  achieved, 
etc?  If  so,  please  elaborate. _ 


6.  General  Comments.  What  do  you  think  should  be  changed  to  improve  future 
reports?  (Indicate  changes  to  organization,  technical  content,  format,  etc.) 


Name 


CURRENT 

ADDRESS 


Organization 


Address 


City,  State,  Zip 

7.  If  indicating  a  Change  of  Address  or  Address  Correction,  please  provide  the 
New  or  Correct  Address  in  Block  6  above  and  the  Old  or  Incorrect  address  below. 


* 


OLD 

ADDRESS 


Name 


Organization 


Address 


City,  State,  Zip 


(Remove  this  sheet  along  the  perforation,  fold  as  indicated,  staple  or  tape 
closed,  and  mail.) 


